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Abstract 

In  many  cases,  optical  tracking  systems  do  not  have  cooperative  beacons  avail¬ 
able.  This  is  particularly  true  for  the  case  involving  tracking  a  laser  illuminated  target 
such  as  a  missile  seeker  head,  where  the  object  of  interest  is  an  extended  source.  Fur¬ 
thermore  the  extended  source  is  often  observed  in  the  presence  of  noise  such  as  shot 
and  speckle  noise  as  well  as  atmospheric  turbulence  which  further  degrades  the  sig¬ 
nal.  This  research  effort  presents  the  evaluation  of  an  existing  algorithm  based  on  the 
maximum- likelihood  technique  for  tilt  estimation  in  the  presence  of  extended  sources 
and  speckle  noise,  with  particular  application  to  the  image  motion  tracking  problem. 
Comparison  is  made  between  the  performance  of  traditional  centroiding  algorithms 
and  the  existing  projection-based  correlation  algorithm  in  simulation.  The  Maximum 
Likelihood  Estimator  using  projection-based  correlation  is  shown  to  offer  improved 
performance  in  the  motion  tracking  problem. 
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Evaluation  of  Performance 
of  a  Maximum  Likelihood  Estimator 
for  Tracking  Purposes  in  the  Presence  of  Speckle  Noise 

I.  Introduction 

This  research  aims  to  evaluate  a  tilt  estimation  algorithm  for  tracking  purposes  that 
performs  well  in  the  presence  of  speckle  noise.  In  particular,  it  addresses  the  type 
of  noise  produced  when  coherent  light  is  reflected  off  an  optically  rough  surface,  as 
distinct  from  the  twinkling  viewed  by  astronomers  which  is  caused  by  atmospheric 
turbulence.  Maximum  Likelihood  Estimation  (MLE)  algorithms  have  been  shown  to 
give  good  performance  compared  to  traditional  centroiding  in  the  presence  of  noise. 
An  MLE  algorithm,  used  in  the  presence  of  speckle  noise,  is  investigated  in  simulation. 
Its  performance  in  terms  of  tilt  error  is  compared  to  that  of  traditional  centroiding. 

1.1  Motivation 

In  the  field  of  Infrared  (IR)  counter  measures,  two  basic  systems  are  employed: 
a  detection  and  tracking  system,  and  a  jamming  system.  On  airborne  platforms  this 
leads  to  several  limitations.  Firstly,  space  and  weight  are  at  a  premium  on  most 
platforms  especially  those  of  the  fighter  type.  Secondly,  with  the  detection  system 
also  providing  tracking  information,  its  ability  to  detect  and  process  new  threats  can 
be  limited.  Thirdly,  a  large  jamming  beam  foot-print  is  often  required,  as  the  location 
in  space  of  the  threat  to  be  jammed  is  not  known  to  any  useful  degree  of  precision. 
This  is  dne  to  the  fact  that  the  exhaust  plume  of  the  threat  is  usually  the  element 
tracked,  not  the  threat  sensor  itself.  These  limitations  can  be  partially  overcome 
if  the  tracking  part  could  be  integrated  with  the  jamming  system.  The  physical 
size  of  the  systems  would  most  likely  be  smaller,  particularly  the  detection  element. 
Once  detection  is  accomplished  and  the  threat  has  been  handed  over  to  the  jamming 
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system,  the  detection  system  can  revert  to  searching  for  new  threats.  Finally,  if  some 
additional  threat  information  is  known,  the  spatial  location  of  the  threat  sensor  to 
be  jammed  can  be  computed  to  a  reasonable  degree  of  accuracy,  allowing  a  smaller 
beam  footprint  and  providing  increased  power  density  in  the  jamming  beam. 

Imaging  wavefront  sensors,  such  as  the  well  known  Shack-Hartmann,  currently 
used  in  the  adaptive  optics  (AO)  field  to  provide  wavefront  subaperture  tilt  estimates, 
can  be  used  to  provide  estimates  of  global  tilt.  This  global  tilt  can  be  used  to  provide 
information  on  a  source’s  position  relative  to  the  optical  axis  of  the  tracking  sensor. 
When  this  information  is  fed  to  a  tracking  system,  the  sensor  can  be  adjusted  such  that 
the  source  is  kept  on  the  optical  axis.  When  this  information  is  also  provided  to  the 
illuminating  source  in  the  case  of  an  IR  jamming  system,  the  efficiency  of  the  jammer 
can  be  much  improved.  Most  sensors  employ  some  form  of  tilt  estimation  algorithm  to 
derive  an  estimate  of  the  global  tilt  parameter.  These  algorithms  generally  perform 
well  in  scenarios  where  the  object  of  interest  (01)  occupies  a  small  portion  of  the 
overall  imaged  scene  and  is  the  single  brightest  feature  in  it,  and  there  is  relatively 
little  to  no  noise  of  any  type.  When  the  first  conditions  is  not  met,  ie.  the  01  occupies 
a  significant  portion  of  the  scene,  the  01  is  termed  non  co-operative.  This  is  the  case 
for  real-world  tracking  scenarios  involving  coherent  illumination,  where  there  is  also 
noise  present  from  the  reflected  light,  the  scene  background,  and  from  within  the 
sensor  itself. 

Few,  if  any,  existing  algorithms  have  been  shown  to  perform  well  in  the  pervi- 
ously  mentioned  scenarios.  An  algorithm  that  performs  well  in  the  presence  of  such 
limiting  factor  would  provide  an  increased  level  of  tracking  accuracy  in  real-world 
situations. 

1.2  Goals 

The  primary  objective  of  this  research  is  to  evaluate  the  performance  of  a 
maximum- likelihood  tilt  estimator  in  the  presence  of  noise,  particularly  speckle.  The 
research  focuses  on  the  performance  of  existing  tilt  estimation  algorithms  in  applica- 
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tion  for  which  they  were  not  originally  intended.  To  satisfy  this  goal,  the  research 
investigates  the  performance  of  an  existing  algorithm  in  ideal  scenarios  as  well  as 
in  non-ideal  scenarios.  The  ideal  scenarios  are  used  to  establish  a  baseline  for  per¬ 
formance  of  each  algorithm.  Non-ideal  effects  are  introduced,  and  the  change  in 
performance  of  each  algorithm  is  examined. 

1.3  Thesis  Outline 

This  thesis  meets  the  goals  through  detailed  analysis  and  simulation.  Chapter  II 
examines  the  theory  behind  the  important  concepts  of  extended  sources  and  speckle 
noise,  including  the  statistics  of  speckle,  upon  which  this  research  is  based.  It  also 
examines  several  of  the  existing  techniques  used  to  estimate  tilt  in  optical  wavefront 
sensing.  Chapter  III  includes  a  detailed  explanation  of  the  structure  of  the  simula¬ 
tions  used  and  reasoning  for  any  assumptions  made.  Chapter  IV  presents  the  results 
of  the  simulation  with  MATLAB®  as  well  as  offering  reasoning  for  the  results.  Chap¬ 
ter  V  provides  concluding  remarks  as  well  as  possible  future  research  opportunities  in 
extending  this  topic. 


3 


II.  Theory  and  Review 


This  work  builds  on  previous  research,  and  moreover,  this  section  reviews  previous 
work  in  the  field  of  tilt  estimation  with  extended  sources.  It  begins  with  the 
basic  theory  of  extended  and  point  sources  and  then  explores  other  important  and 
relevant  areas  that  are  the  foundations  of  this  research. 

2.1  Point  and  Extended  Sources 

While  point  sources  seem  at  first  glance  to  be  the  simplest  beacons  to  use  in  sim¬ 
ulation,  this  is  in  fact  not  entirely  true.  If  one  were  to  accurately  model  a  point  source, 
it  would  have  to  be  represented  as  a  Dirac  delta  function  in  space.  It  has  infinite  band¬ 
width,  which  is  clearly  not  possible  in  simulation  because  of  a  discretely  sampled  grid 
and  impossible  in  a  real- world  situation.  Firstly  we  attempt  to  define  the  physical 
differences  between  point  and  extended  sources.  When  trying  to  characterize  point 
and  extended  sources,  it  is  often  helpful  to  start  with  the  geometrical  perspective.  A 
true  point  source  at  a  finite  distance  emits  divergent  light  in  all  directions,  while  at 
very  large  distances  light  from  a  point  source  can  be  considered  collimated  [12],  An 
extended  source  produces  overlapping  pencils  of  rays  called  beams  from  each  object 
point.  Figures  2.1  and  2.2  illustrate  the  spreading  effect  of  an  extended  source. 

A  point  source  can  be  thought  of  in  a  number  of  ways.  One  is  to  simply  refer 
to  it  as  resolved  or  unresolved,  in  that  it  either  occupies  more  than  one  pixel  in  the 
image  sensor  or  not.  However,  this  may  not  capture  the  physics  of  any  propagation 
involved.  Another  is  an  idealized  source  whose  dimensions  are  very  small  compared 
to  the  viewing  distance,  i.e.  the  product  of  the  lateral  source  dimensions  should  be 
very  much  smaller  than  the  distance  to  the  detector  or  observation  screen  squared 
and  the  included  solid  angle  is  very  small.  Hence  we  can  define  a  point  source  as 
shown  Figure  2.3. 

Radiant  energy  emitted  by  a  point  source  is  referred  to  as  isotropic,  i.e.  it 
radiates  equally  in  all  directions.  The  surface  area  through  which  the  radiation  passes 
is  47rr2  at  a  distance  r  from  the  source.  If  a  detector  of  area  A  is  placed  a  distance  r 
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Figure  2.2:  Light  rays  from  a  distant  extended  object. 
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Figure  2.3:  Defining  a  point  source. 


from  the  source,  the  detector  subtends  a  solid  angle  D  =  A/r2.  The  radiant  intensity 
can  be  expressed  by 


Ie  = 


4*e 

TT’ 


(2.1) 


where  0e  is  the  radiant  flux,  defined  as  the  measure  of  the  total  power  of  electro¬ 
magnetic  radiation  (including  infrared,  ultraviolet,  and  visible  light)  landing  on  a 
particular  surface.  The  measured  intensity  remains  constant  as  the  detector  moves 
away  from  the  source.  However,  we  can  rarely  neglect  diffraction  effects  and  the  mea¬ 
sured  intensity  takes  the  form  of  Ig  =  Ig( 0)  cos3#  from  Ref.  [17],  illustrated  in  Figure 
2.4. 


Extended  sources  require  special  treatment  to  properly  measure  corresponding 
flux  and  surface  brightness.  This  arises  from  the  complex  scattering  of  incident  light 
reaching  the  sensor  focal  plane  array.  The  images  of  extended  sources  have  more 
extended  wings  than  those  of  point  sources.  Photons  that  would  normally  be  scattered 
out  of  the  point  spread  function  (PSF)  aperture  used  to  measure  a  point  source 
are  instead  captured  when  an  extended  source  is  present.  Extended  sources  can  be 
considered  collection  of  point  sources  as  suggested  in  Ref.  [19].  For  an  extended  source, 
the  point  source  model  needs  to  be  expanded  to  account  for  the  area  of  the  source 
and  the  fact  that  the  radiant  intensity  may  vary  across  this  area.  If  we  consider  the 
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Figure  2.4:  Variation  in  irradiance  from  a  point  source  on  a  detector  at  distance  r 


extended  source  to  comprise  equally  spaced  point  sources  as  shown  in  Figure  2.5(a), 
then  the  irradiance  per  unit  area  Le  can  be  expressed  as 


Le  = 


4*e 

An1 


(2.2) 


where  <pe  is  the  radiant  flux  of  the  source,  A  is  the  area  of  the  source  and  is  the 
solid  angle  subtended  by  the  detector  at  the  source.  If  the  same  extended  source  is 
now  viewed  at  an  angle  6  as  shown  in  Figure  2.4,  then  the  spacing  between  point 
sources  is  shortened  in  one  direction  and  the  area  A  is  reduced  by  a  factor  cos  6  as  in 
Figure  2.5(b).  Thus  A  is  replaced  by  dAcosd  in  Eq.  (2.2)  leads  to 

LM  =  (2.3) 

cost/ 

where  Le( 0)  is  the  radiance  at  6  =  0.  Unfortunately  with  Eq.  (2.3)  we  are  about 
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Figure  2.5:  Extended  source  considered  as  a  collection  of  point  sources  when 

viewed(a)  at  normal  to  the  source  and  (b)  from  an  angle  9  to  normal. 


to  violate  a  few  physical  laws  in  that  it  suggests  that  as  6  decreases  the  radiance 
increases,  eventually  to  infinity  at  9  =  n/2.  One  way  to  avoid  this  is  to  choose  an 
extended  source  whose  radiance  is  independent  of  viewing  angle.  This  is  possible  if 
the  radiant  intensity  falls  off  as  cos  9  which  compensates  for  the  effect  in  Eq.  (2.3), 
which  now  becomes 

le(9)  =  cos9j^-  (2-4) 

The  extended  source  can  be  thought  of  as  a  lambertian  source  as  shown  in 
Figure  2.6  whose  intensity  varies  as  the  cosine  between  the  normal  and  the  direction 
of  propagation.  It  should  be  noted  here  that  most  practical  sources  do  not  strictly 
meet  this  requirement,  though  many  are  close  enough  to  make  the  approximation 
valid.  Now  to  calculate  the  irradiance  at  a  distance  r  from  the  extended  source,  we 
note  the  irradiance  produced  by  an  element  dA  of  the  source  is  given  by  Eq.  (2.2) 
which  when  substituted  into  Eq.  (2.1)  and  after  a  bit  of  manipulation,  gives  the 
irradiance 


Len. 


E, 


(2.5) 


Figure  2.6:  Lambertian  source  where  the  intensity  reduces  as  the  cosine  of  the  angle 
between  the  normal  and  the  direction  of  propagation. 


The  question  then  is,  when  is  it  no  longer  accurate  to  approximate  a  physical  source 
by  a  point  source?  If  we  take  the  circular  source  irradiance  as 


Ep  = 


(2.6) 


where  a  is  the  source  radius  and  r  is  the  detector  distance,  the  extended  disk  source 
irradiance  is 


_  n  a2L 

o  9  ■ 

a  +  r 


(2.7) 


Now  if  we  desire  to  have  less  than  1%  error  by  approximating  the  disk  as  a  point  then 


EP 


EP 


<  0.01  a  <  0.01  x  r 


(2.8) 


This  is  the  basis  for  the  “5  times”  rule  of  thumb,  which  is  the  observation  distance 
should  be  at  least  five  times  the  largest  source  dimension.  The  reader  can  compare  this 
to  the  definition  presented  earlier  in  this  section.  The  physical  differences  between 


extended  and  point  sources  reduce  wavefront  sensing  accuracy  by  various  degrees 
depending  on  the  sensor  used. 

2.2  Atmospheric  Turbulence 

Random  variations  in  the  temperature  and  pressure  of  the  earth’s  atmosphere 
alter  the  the  refractive  index  of  the  air,  both  spatially  and  temporally.  Optical  waves 
propagating  through  these  variations  are  distorted.  This  distortion  is  known  as  tur¬ 
bulence.  Turbulence  affects  all  optical  systems  which  propagate  light  through  long 
atmospheric  paths.  There  have  been  many  theories  and  much  work  done  on  char¬ 
acterizing  the  effects  of  turbulence  on  optical  propagation.  A  statistical  analysis  is 
necessary  as  it  is  not  feasible  to  describe  exactly  the  refractive  index  at  all  points  in 
space  and  time.  The  most  widely  accepted  theory,  due  to  its  consistency  with  obser¬ 
vation,  is  that  put  forward  by  A.N.  Kolmogorov  [11],  His  theory  centers  on  randomly 
distributed  pockets  of  air,  called  eddies,  of  varying  sizes  and  temperatures  causing  the 
random  variations  in  the  refractive  index  of  the  atmosphere.  From  this,  a  refractive 
index  profile  of  the  section  of  atmosphere  in  question  can  be  developed. 

Kolmogorov  suggested  that  turbulent  flow,  governed  by  the  Navier-Stokes  equa¬ 
tions,  could  be  described  by  the  transfer  of  kinetic  energy  from  large  eddies  into  smaller 
eddies.  The  average  size  of  these  eddies  being  designated  the  outer  scale  L0  for  the 
large  eddies  and  the  inner  scale  /0  for  the  smaller  eddies.  L0  can  range  from  the  height 
above  ground  at  lower  altitudes  up  to  hundreds  of  meters  at  higher  altitudes.  /0  ranges 
from  a  few  millimeters  at  low  altitudes  to  centimeters  higher  up.  The  range  of  sizes 
between  inner  and  outer  scales  is  known  as  the  inertial  subrange,  and  Kolmogorov 
assumed  that  eddies  within  this  range  are  statistically  homogeneous  and  isotropic.  It 
is,  however,  more  accurate  to  say  that  properties  such  as  refractive  index  and  wind 
velocity  have  stationary  increments.  Kolmogorov  determined  that  the  average  speed 
v  of  turbulent  eddies  is  related  to  their  scale  size  r  via 


Kolmogorov’s  analysis  of  potential  temperature  (linearly  related  to  regular  temper¬ 
ature  T)  led  to  an  expression  for  the  refractive  index  at  a  point  in  space  r  given 
by 

n(r)  =//n(r)  +  ni(r),  (2.10) 

where  fin(r)  is  the  slowly  varying  mean  of  the  refractive  index  value,  and  7ii(r)  is 
the  deviation  of  the  index  from  its  mean  value.  This  creates  a  zero-mean  random 
process  ni(r).  At  optical  wavelengths  considered  here,  the  refractive  index  of  air  can 
be  approximated  by 


P(y) 

n(r)  =  1  +  7.99  x  10~5-f4  for  A  =  0.5 jum, 

T(r) 


(2.11) 


where  A  is  the  optical  wavelength,  P  is  the  pressure  in  millibars  and  T  is  temperature 
in  Kelvin.  Then  assuming  each  eddy  has  relatively  uniform  pressure,  the  variation  in 
the  refractive  index  is  given  by 


dn 


7.99  x  1(T5 


dB 

rp'2  ' 


(2.12) 


Variation  in  refractive  index  is  directly  proportional  to  the  variation  in  potential 
temperature  and  as  such,  the  refractive  index  structure  function  follows  a  similar 
power  law  to  Eq.  (2.9): 


Dn( r)  =  C'?y2/:!  for  l0  <  r  <  L0 


(2.13) 


where  C'2  is  the  refractive-index  structure  function  parameter,  in  m  2/3. 

It  is  often  necessary  to  have  a  spectral  description  of  the  fluctuation  in  refrac¬ 
tive  index.  The  Kolmogorov  power  spectral  density  $n(/c)  can  be  computed  from 
Eq.  (2.13)  and  is  given  by 

<f>n(«)  =  0.033C^tt-11/3  for  y—  <  k  <  (2-14) 

A0  <0 
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where  k  =  2tt  (^fx i  +  is  the  angular  spatial  frequency  in  rad/m. 

Through  the  application  of  Rytov  theory  (Maxwell’s  equations  and  perturbation 
theory),  the  field  can  be  written  as 

U(r)  =  C/o(r)exp['0(r)],  (2.15) 

where  ?7o(r)  is  the  vacuum  solution  (n i  =  0)  of  Maxwell’s  equations  and  if( r)  is  the 
complex  phase  perturbation.  Successive  perturbations  of  the  form 

ip(r)  =  -0i(r)  +  ip2{r)  +  ■  ■  ■  (2.16) 

can  be  used  to  compute  statistical  moments  of  ^  which  in  turn  are  used  to  compute 
statistical  moments  of  the  field.  For  example,  turbulent  phase  screens  such  as  those 
used  in  simulation  in  Section  3.2.4  are  realizations  of  ^(r). 

Rytov  theory  yields  many  parameters  that  can  be  used  to  characterize  optical 
impact  of  turbulence.  Two  of  these  are  the  coherence  diameter  ro  and  the  isoplanatic 
angle  6q.  both  of  these  parameters  are  computed  from  the  integrated  moments  of  the 
structure  parameter  Cf.  Cf  is  the  measure  of  the  local  turbulence  strength  at  a  point 
in  space,  and  is  a  function  the  propagation  distance  Az.  Functions  of  C'2,  known  as 
structure  functions ,  describe  the  local  turbulence  along  a  particular  optical  path.  The 
coherence  diameter  r0,  also  known  as  the  Fried  parameter,  is  given  approximately  by 

/  \  5/3 

A/,( r)  =  6.88  (  -)  .  (2.17) 

Values  of  r o  are  typically  5-10  cm  at  visible  wavelengths.  The  isoplanatic  angle  do  is 
defined  as  the  angle  between  two  point  sources  for  which  the  mean  square  phase  differs 
by  1  rad2.  It  may  also  be  considered  as  the  largest  field  angle  over  which  the  optical 
path  through  the  turbulence  does  not  vary  significantly  from  the  on-axis  optical  path 
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through  the  turbulence.  Values  for  90  at  visible  wavelengths  are  typically  5  —  lOyurad, 
looking  directly  overhead. 

If  the  turbulence  encountered  by  the  beam  is  strong,  then  the  beam  from  an 
extended  source  wanders  over  time.  If  a  coherent  source  such  as  a  laser  is  used 
to  illuminate  a  rough  target,  the  beam  also  exhibits  random  effects  which  include 
strong  intensity  and  phase  variations.  These  variations  create  challenging  scenarios 
for  wavefront  sensing.  Light  rays  arriving  at  a  sensor  aperture,  such  as  a  Shack- 
Hartmann  wavefront  sensor,  from  different  points  of  an  extended  source  have  traveled 
significantly  different  atmospheric  paths.  Light  corrupted  by  aberrations  obtained 
from  these  different  paths  arrives  superimposed  at  the  aperture  where  it  is  used  in 
the  sensing  system.  Considering  the  light  is  already  superimposed  and  has  most 
probably  arrived  from  many  different  directions,  often  separated  by  more  than  90,  the 
isoplanatic  angle,  given  by  , 


90  =  (L09fc2C^L8/3)5/3,  (2.18) 

for  a  constant  C2  path,  or  from  points  separated  by  distances  exceeding  ro,  the  co¬ 
herence  diameter,  given  by 

Az 

r0,sw  =  0.423/c2  y  C^(z)(-^)5/3dz 
o 

conventional  wave-front  sensing  processes  do  not  estimate  the  the  same  wave-front 
error  that  would  be  due  to  a  point  source.  In  many  cases  the  error  in  estimation  is 
significant  and  severely  degrades  the  sensor’s  ability  to  correctly  estimate  the  wave- 
front.  Nominally,  when  an  extended  source  is  present,  the  intensity  pattern  in  the 
pupil  of  the  observation  system  should  be  a  lower-contrast,  blurred  version  of  the 
intensity  pattern  expected  of  a  point  source. 

Often  a  Hartmann-type  sensor  is  used  for  wavefront  sensing.  The  presence  of 
an  extended  source  in  this  case  also  creates  some  interesting  effects.  If  the  complex 


(2.19) 
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Figure  2.7:  Typical  Hartmann  lenslet  array. 

transmittance  of  a  Hartmann  lenslet  array  such  as  the  one  shown  in  Figure  2.7  is 
given  by 


tH{xP)  =  J^exp [-j^y(xp  ~  xs)2]rect  fXp  ^  ,  (2.20) 

s 

where  S  is  the  number  of  sub-apertures,  xs  is  the  center  of  the  sth  sub-aperture,  fi 
is  the  focal  length,  d  is  the  side  dimension  of  the  lenslet  and  k  is  the  wave  number. 
When  a  field  due  to  an  extended  source  falls  on  a  Hartmann  array,  due  to  the  nature 
of  the  light  from  the  extended  source,  some  held  segments  may  overlap  onto  adjacent 
sub- apertures.  The  resultant  intensity  is  derived  in  Ref.  [21]  and  is  shown  here  as 


Ih(xHi  t ) 


dxpdx'pdxphpiF^xp^  XH)h*NF(x'P,  XH)tH(xp)t*H(x'P ) 


'xp 


’  XT 


x  | UT(xT,t)\2hA(xT,xp,t)h*A(xT,x'p,t),  (2.21) 
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Conventional  Image,  Point  Source  Beacon  Case 


Conventional  Image,  Random  Extended  Beacon 
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Figure  2.8:  (a)  Conventional  image  point  source  case,  (b)  Conventional  image 

extended  beacon  case.  [21] 


where  hNp(xp,  %h)  is  the  impulse  response  for  the  propagation  from  the  lenslct  array 
to  the  sensor  detector  plane,  Xp  is  a  coordinate  in  the  Hartman  detector  plane  and 
h*A(xT,  x'P,t)  is  the  impulse  response  of  the  atmosphere  [21].  The  integrals  over  xp 
and  x'p  are  propagations  which  move  the  held  back  from  the  lenslet  array  to  the 
detector  plane  and  their  product  becomes  the  detector  intensity.  The  integral  over  xp 
would,  under  isoplanatic  conditions,  be  a  convolution  which  describes  the  propagation 
of  the  held  from  the  beacon  to  the  aperture.  However  when  considering  an  extended 
beacon  and  possibly  anisoplanatic  conditions  associated  with  it,  the  Hartmann  sensor 
has  a  significantly  different  result  than  for  just  a  point  source.  This  is  illustrated  in 
Figure  2.8  with  images  from  Ref.  [21]  where  image  (a)  clearly  has  a  different  centroid 
location  from  image  (b). 


2.3  Speckle  Noise 

2.3.1  Speckle  Phenomena.  When  complicated  objects  are  illuminated  by 
highly  coherent  light  of  the  type  produced  by  a  laser,  an  important  type  of  image 
defect  is  seen.  In  particular,  whenever  the  object  is  rough  on  the  scale  of  an  optical 
wavelength,  the  image  has  a  grainy  appearance.  The  contrast  in  the  image  is  very 
pronounced,  with  large  numbers  of  bright  and  dark  spots.  These  spots  apparently 
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Figure  2.9:  Typical  Speckle  Noise  Pattern. 

have  relationship  to  the  scattering  properties  of  the  surface  of  the  scattering  object. 
This  image  defect  is  called  speckle.  The  effect  is  not  limited  to  reflections  from 
objects  being  also  seen  in  images  of  transparent  objects  illuminated  by  coherent  light 
through  a  diffuser.  A  typical  speckle  pattern  is  shown  in  Figure  2.9.  Detailed  analysis 
of  speckle  began  in  the  1960’s.  However,  studies  had  actually  been  carried  out  far 
earlier  (although  not  known  as  speckle  at  the  time),  by  Verdet  (1865)  and  Lord 
Rayleigh  (1880)  in  their  work  on  Fraunhofer  rings.  Von  Laue  [14-16]  derived  many 
of  the  basic  speckle  properties  in  the  study  of  light  scattered  from  a  large  number  of 
particles. 

The  vast  majority  of  surfaces  are  extremely  rough  on  the  scale  of  an  optical 
wavelength  [9]  and  under  illumination  by  monochromatic  light,  the  wave  reflected 
form  such  a  surface  consists  of  contributions  from  many  scattering  points  with  random 
phases.  The  image  formed  at  a  given  point  in  an  observation  plane  consists  of  a 
multitude  of  amplitude  spread  functions,  each  arising  from  a  separate  point  on  the 
rough  scattering  surface  as  shown  in  Figure  2.10.  As  a  result,  the  contributing  spread 
functions  have  widely  varying  phase  and  when  added  together,  produce  a  highly 
complex  interference  pattern. 

The  same  arguments  apply  to  transmission  objects  illuminated  via  a  diffuser. 
The  diffuser  causes  the  exiting  wavefront  to  have  a  highly  complex  and  corrugated 
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Object 

Figure  2.10:  Speckle  formation  in  the  image  of  a  rough  object,  (redrawn  from  [9]) 

structure.  In  the  image  of  this  type  of  object,  we  can  again  see  the  large  intensity 
fluctuations  caused  by  overlapping,  dephased  amplitude  spread  functions  [9]. 

2.3.2  Speckle  Intensity  Statistics.  Our  knowledge  of  the  exact  wavelength- 
scale,  structure  of  the  wavefront  leaving  the  surface  is  often  extremely  limited,  so  it 
is  beneficial  to  think  in  terms  of  the  statistical  properties  of  speckle  patterns.  The 
statistics  are  defined  over  an  ensemble  of  objects,  all  with  the  same  properties  but 
differing  in  the  detail.  If  a  detector  is  placed  in  the  image  plane  at  a  precisely  known 
position,  the  measured  intensity  cannot  be  predicted  exactly.  We  therefore  attempt 
to  predict  the  intensity  by  using  the  statistical  properties  of  the  intensity  over  an 
ensemble  of  rough  surfaces.  A  most  important  statistical  property  of  speckle  is  the 
probability  density  function  (PDF)  of  the  observed  intensity  /  at  a  point  in  the  image. 
Effectively  we  are  asking,  how  likely  is  it  that  we  will  see  a  spot  of  given  intensity  at 
that  point?  It  is  analogous  to  the  well-known  random  walk  [23,24],  If  the  phases  of 
the  individual  scattered  contributions  are  uniformly  distributed  over  (— 7T,  7 r),  i.e.  the 
object  is  rough  on  the  scale  of  a  wavelength,  then  the  field  associated  with  any  single 
linear  polarization  component  of  the  the  image  is  a  circular  complex  Gaussian  random 
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I  =  5  W/m2 

Figure  2.11:  Negative  Exponential  distribution  with  mean  7=5. 
variable  (RV),  and  its  intensity  statistics  have  a  negative  exponential  distribution, 

{7  exp  (~)  for  I  >  0, 

'  1 

0  elsewhere 

where  7  is  the  mean  intensity  associated  with  that  polarized  component.  Figure  2.11 
shows  an  example  of  such  a  negative  exponential  distribution.  If  the  scattered  wave  is 
partially  polarized,  it  can  be  shown  that  the  density  function  for  /  is  the  difference  of 
two  negative  exponential  functions  [9].  The  negative  exponential  distribution  nature 
of  the  speckle  intensity  implies  the  fluctuations  around  the  mean  are  very  pronounced. 
If  we  define  the  contrast  of  a  speckle  pattern  as  the  ratio  of  standard  deviation  to 
mean,  which  for  the  polarized  case,  C  —  ai/I  —  1.  It  is  because  of  this  very  high 
contrast  that  speckle  is  extremely  disturbing  to  the  human  eye. 

It  is  important  to  mention  that  the  distribution  of  mean  intensity  I(x,y)  in 
the  image  of  a  coherently  illuminated  object  is  identical  to  the  image  intensity  that 
would  be  observed  if  the  object  were  illuminated  with  spatially  incoherent  light  of 
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p, 

p2 

Figure  2.12:  Typical  Young’s  Experimental  Setup  (redrawn  from  [9]) 

the  same  spectral  power  density  [9].  Incoherent  illumination  may  be  regarded  as  a 
rapid  time  sequence  of  spatially  coherent  wavefronts.  Thus  the  time-integrated  image 
intensity  observed  from  spatially  incoherent  illumination  is  identical  to  the  ensemble 
average  intensity  I(x,y)  (assuming  identical  bandwidth).  Therefore,  methods  used 
for  incoherent  intensity  distribution  may  also  be  used  for  predicting  the  mean  speckle 
intensity  distribution  with  coherent  illumination  of  a  rough  object. 

When  optically  rough  objects  are  illuminated  by  monochromatic  light,  the  re¬ 
flected  wave  cannot  be  treated  as  ergodic  [9]  as  the  time  and  ensemble  averages  are  not 
equal.  This  can  be  demonstrated  by  considering  two  different  Young’s  experiments  of 
the  type  shown  in  Figure  2.12. 

First,  let  light  reflected  from  a  rough  surface  fall  on  a  mask  containing  two  pin¬ 
holes.  The  fringe  formed  is  observed  on  a  distant  screen.  As  the  light  is  monochro¬ 
matic,  it  is  also  spatially  coherent  [9]  and  the  fringe  has  visibility  v  given  by 

2  Vhh 

v  =  t  T7  ’ 

where  1 1  and  1-2  are  the  intensities  of  light  falling  on  pinholes  one  and  two,  respec¬ 
tively.  Since  the  distribution  does  not  change  with  time  and  the  coherence  of  the 
light  has  not  been  reduced  by  the  amplitude  and  phase  distribution  imparted  onto 
the  wave  by  a  rough  surface,  the  modulus  of  the  complex  coherence  factor  /i 1 2 1,:  as 
defined  in  Table  2.1,  must  be  one  (time  averaged  definition  of  coherence).  In  the 
second  example,  objects  with  different  surface  profiles  are  successively  placed  in  the 
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Table  2.1:  Names  and  definitions  of  Various  measures  of  Coherence  used  in  reference  to  Figure  2.12  (reproduced  from  [9]) 


Symbol 

Definition 

Name 

Temporal  or  Spatial  Coherence 

fuM 

[Note:rn(0)  =  /(Pi)] 

Self  coherence  function 

Temporal 

7n(T) 

I’n  C7"  1 
rn(o) 

Complex  degree  of  (self)  Coherence 

Temporal 

Tl2(T) 

{u{P1,t  =  TU*(P2,t)) 

Mutual  coherence  function 

Spatial  and  temporal 

712  (t) 

l'i2  (r) 

frn(o)r22(o)i1/2 

Complex  degree  of  coherence 

Spatial  and  temporal 

J 12 

(m(Pi,i  =  rM*(Pi,t))  =  r12(o) 

Mutual  intensity 

Spatial  quasimonochromatic 

H 12 

ir7>  -  •*««» 

Complex  coherence  factor 

Spatial  quasimonochromatic 

illuminating  beam,  and  all  of  the  generated  fringes  are  time  integrated.  Any  one  of 
these  component  fringes  has  a  visibility  corresponding  to  \]JLyi\  —  1,  but  the  super¬ 
position  of  the  successive  fringes  does  not  because  the  phase  from  fringe  to  fringe  is 
not  constant  as  each  amounts  to  a  new  realization.  Therefore  the  ensemble  averaged 
fringe  produces  a  \]ivi\  that  is  not  equal  to  1. 

The  wave  equation  governing  the  propagation  of  light  remains  the  same  for 
time-  or  ensemble-averaged  properties  of  light.  As  such  the  laws  governing  the  propa¬ 
gation  of  coherence  functions  for  time-  and  ensemble-averaged  qualities  are  identical. 
Therefore  the  mutual  intensity,  as  defined  in  Table  2.1,  of  light  reflected  from  a  rough 
surface  and  observed  very  close  to  that  surface  is  the  same  as  that  observed  from  an 
incoherent  source.  Over  an  ensemble  of  ideally  rough  surfaces  there  is  little  relation¬ 
ship  between  the  phases  of  scattered  light  from  two  closely  spaced  (<  A)  elements  on 
the  surface  represented  by  a  delta  function  mutual  intensity, 


i;  62,^72)  =  kI(Zi,v i)<K6  -6,  0  -V2), 


where  /  is  the  ensemble  average  intensity  distribution  and  k  is  a  constant.  The  mutual 
intensity  observed  at  a  distance  z  from  the  source  can  be  computed  using  the  Van 
Cittert-Zernike  theorem  and  is  given  by 


J(x1,y1;x2,y2) 


ne 


+OO 


(A^) 


27T 

/(£,  rj)  exp  \  j  —  [(Axt  +  Ayr])} 


where  ^  =  j-\{x\  +  y|)  —  {x\  +  yf)],  (£,77)  is  a  point  in  the  source  plane,  (x,y)  is 
a  point  in  the  observation  region,  Ax  =  x-i  —  aq,  Ay  =  y2  —  yi  ,  A  is  the  mean 
wavelength  of  the  source  and  /(£,y)  is  the  ensemble  averaged  intensity  distribution 
across  the  scattering  spot  on  the  rough  surface.  Given  the  geometry  as  in  Figure  2.10, 
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the  mutual  intensity  in  the  image  is  given  by 


J(ui,v1;u2,v2) 


dxdy , 


—  OO 

where  P  is  the  complex  pupil  function  of  the  lens,  n  is  a  constant  and  A u  =  w2  — 

Ml,  An  =  1>2  —  Ml- 


2.3.3  Speckle  Photon  Count  Statistics.  While  intensity  distribution  statistics 
are  both  helpful  and  valuable  aids  in  describing  speckle,  photon  count  statistics  are  at 
the  heart  of  detection  and  motion  estimation  for  speckle  images.  We  consider  a  single- 
mode  laser  whose  light  falls  on  a  detector  and  we  wish  to  determine  the  distribution  of 
the  number  of  events  in  any  r-second  interval.  Assuming  constant  incident  intensity, 
the  integrated  intensity  on  a  pixel  with  area  A  is  given  by 

W  =  IqAt, 


with  probability  density  of  the  form 

Pw(W)  =  8(W  -  I0At). 

Substituting  into  Mandel’s  formula  [18],  in  which  the  unconditional  probability  of 
observing  K  photo  events  can  be  expressed  as 


P(K)  —  /  P(K\W)pw  (W)dW 

Jo 

f°°  ( aW)K 


(2.22) 


K\ 


exp~aW  pw  (W)  dW, 


where  a  is  a  proportionality  constant  equal  to  the  ratio  of  quantum  efficiency  to 
energy  rj/hv  and  pw  (W)  is  the  probability  density  function  of  the  integrated  intensity. 
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Performing  the  integration  yields 


Pk{K) 


(aI0Ar)K  InAr 
K\ 


(2.23) 


Since  the  mean  number  of  photon  events  K  =  cd0Ar,  we  can  rewrite  Eq.  (2.23)  as 


Pk{K) 


(Kfe-K 

K\ 


From  this,  the  mean  in  terms  of  the  variance  is  a2K  =  K .  With  care,  this  ideal  model 
can  be  closely  approximated  in  practice.  There  are  two  cases  of  photon  count  statistics 
for  polarized  thermal  radiation: 

1.  Counting  time  shorter  than  coherence  time,  and 

2.  Counting  for  arbitrary  time. 

we  begin  with  the  first  case,  it  gives  a  good  approximation  for  the  tracking  scenario. 
Since  the  counting  time  is  extremely  short,  the  incident  intensity  can  be  assumed  to 
be  constant  over  the  entire  counting  interval.  Therefore,  the  integrated  intensity  is 
equal  to  the  product  of  the  intensity,  the  counting  time,  and  the  detector  area: 


W  =  I  {t)  At. 


The  value  of  the  intensity  within  that  interval  is  random  and  obeys  negative  expo¬ 
nential  statistics  as  discussed  in  Section  2.3.2.  It  follows  that  the  same  should  be  true 
of  the  integrated  intensity,  such  that 

Pw{-W)  =  =  exp  (-^-) .  w  a  o. 


We  can  now  find  the  photon  count  statistics  by  substituting  into  Mandel’s  equation 
and  integrating  to  give 


P(K) 


1  /  aW  \ 

1  +  aW  \1  +  aW  ) 
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Substituting  in  K  =  aW  gives 


P(K)  =  —^=(— .  (2.24) 

V  '  1  +  K  V 1  +  K  J  K  ' 

The  distribution  represented  by  Eq.  (2.24)  is  the  Bose-Einstein  distribution  which  has 
variance  equal  to  a2K  =  K  +  (A')2.  Note  that  the  first  term,  K ,  represents  the  Poisson 
nature  of  the  integration  of  light  and  matter  and  the  second,  (A')2,  represents  the 
fluctuations  of  integrated  intensity  which  are  significant  if  K  3>  1.  If  we  examine  the 
signal-to-noise  ratio  given  by 

S  _  I  K 
N~  V  i+7T 

it  can  be  seen  that  the  S/N  approaches  1  as  K  — >  oo.  This  indicates  that  the  count 
variation,  just  as  with  the  intensity  variation,  is  substantial.  Probability  masses 
associated  with  Poisson  and  Bose-Einstein  distributions  are  shown  in  Figures  2.13 
and  2.14,  respectively.  Comparison  of  these  two  figures  shows  that  when  the  mean 
number  of  counts  is  greater  than  1,  the  spread  of  the  Bose-Einstein  distribution  is 
greater  than  that  of  a  Poisson  distribution  and  consequently  fluctuations  in  photon 
count  for  Bose-Einstein  is  greater  than  for  Poisson.  Additionally,  when  the  number 
of  counts  K  is  -C  1,  the  differences  between  the  two  distributions  is  small,  and  it  can 
be  shown  [9]  that  only  one  and  zero  events  have  significant  probability  so  that  the 
two  distributions  become  asymptotically  the  same. 

Alternatively,  since  this  research  is  giving  consideration  to  shot  noise  as  well  as 
speckle  noise,  the  negative  binomial  distribution  case  shown  by  Goodman  in  Ref.  [9] 
is  also  examined.  This  case  includes  the  effects  of  both  the  random  arrival  nature  of 
photons  and  the  negative  exponential  distribution  of  speckle  noise.  The  difference  in 
this  case  is  that  the  counting  interval  r,  previously  assumed  to  be  much  shorter  than 
the  coherence  time  of  the  incident  light,  is  now  an  arbitrary  interval  which  maybe 
longer  than  the  coherence  time.  Assuming  that  the  wave  incident  on  the  sensor  has 
a  coherence  area  that  is  much  larger  than  the  area  of  the  sensor,  the  approximate 
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0.2 


K  =  5 

Figure  2.13:  Poisson  distribution  with  mean  K  =  5 


Figure  2.14:  Bose-Einstein  distribution  with  mean  K  =  5 
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solution  for  the  probability  density  Pw(W)  is  a  gamma  probability  density  [9],  given 
by 


Pw(W) 


M  exp(— _M  =  ) 

I XM) 


for  W  >  0, 


elsewhere 


where  A4  is  the  number  of  degrees  of  freedom  of  the  intensity  included  in  the  mea¬ 
surement  interval  (see  Ref.  [9]  for  more  detail).  Knowing  the  approximate  form  of  the 
probability  density  function  of  the  integrated  intensity,  the  probability  density  func¬ 
tion  of  the  number  of  photo  counts  in  the  arbitrary  time  interval  can  be  calculated. 
Again  using  Mandel’s  formula  and  integrating  gives 


P(K) 


r(K  +  M) 
r(K  +  i)r(M) 


_  M 

-K 

K 

1  + 

1  +  — 

A 

M 

-M 


(2.25) 


where  K  =  aW.  This  is  known  as  the  negative  binomial  distribution  and  is  a  good 
approximation  [9]  to  the  photo  count  distribution  considered  in  this  research. 


2-4  Tilt  Estimation 

Wave  fronts,  also  called  phase  fronts  can  be  described  as  a  line  along  which  all 
points  have  the  same  phase,  i.e.  a  surface  of  constant  optical  path  length  (OPL).  It 
would  be  fantastic  if  we  could  directly  measure  this  wavefront,  however,  at  the  visible 
and  infrared  frequencies  concerned,  the  phase  of  the  light  does  not  interact  with  the 
medium  through  which  it  travels  in  a  manner  which  we  can  observe.  Just  as  our  eyes 
respond  to  changes  in  intensity,  most  detectors  also  respond  to  the  intensity  of  the 
incident  light.  We  need  a  way  of  using  this  fact  to  derive  the  phase  of  the  wavefront 
from  the  observed  change  in  intensity.  There  are  many  types  of  aberrations  resulting 
from  changes  imposed  on  the  wavefront.  We  are  only  concerned  here  with  tilt  and 
the  measurement  of  it  for  tracking  purposes.  Tilt  can  be  described  as  the  deviation 
of  the  incident  wavefront  from  a  reference. 

Figure  2.15  shows  how  a  beam  with  a  tilted  wavefront,  passing  through  an 
aperture  follows  a  set  geometry.  If  the  beam  is  focussed  to  a  spot,  and  the  beam 
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Figure  2.15:  Wavefront  Tilt  Geometry  -  redrawn  from  [26] 

has  no  tilt,  then  the  spot  would  be  focused  on  the  optical  axis.  If  the  beam  has 
tilt  then  the  spot  is  focussed  off  axis.  The  spot  position  is  shifted  off  axis  by  a 
distance  proportional  to  the  tilt  angle.  In  actual  fact,  it  is  the  weighted  centroid  of 
the  focussed  spot  that  is  shifted.  If  a  detector  which  responds  to  the  position  of  the 
focussed  spot  is  placed  in  the  focal  plane  we  have  a  way  of  measuring  the  tilt.  Tilt 
sensors  measure  the  OPL  difference  either  directly  in  the  form  of  an  interferogram  or 
indirectly  through  the  differential  wavefront  as  a  function  of  pupil  coordinates.  Tilt 
sensors  convert  angular  wavefront  errors  into  intensity  variations  that  can  be  sensed 
by  a  photodetector  and  converted  to  a  wavefront  tilt  measurement  by  other  means. 
The  traditional  tilt  estimation  and  motion  tracking  problem  is  well- understood  and 
adequately  addressed  when  the  object  of  interest  is  a  point  source  or  a  source  whose 
dimensions  and  distance  from  the  observing  aperture  allow  it  to  be  safely  treated 
as  a  point  source.  However,  when  the  dimensions  of  the  object  of  interest  preclude 
it  from  being  safely  modeled  as  point  source,  as  in  many  real-world  scenarios,  the 
problem  is  not  so  well-addressed.  The  uncooperative  or  extended  source  presents 
many  challenges  to  the  motion  tracking  and  tilt  estimation  problem.  The  traditional 
centroiding  method  employed  in  most  Shack- Hartmann  type  sensors  relies  on  points 
of  high  contrast  within  the  image  to  act  as  beacons  in  estimating  the  wavefront 
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tilt.  When  considering  extended  beacons  or  scenes  which  are  large  compared  to  the 
aperture  diameter,  there  are  often  few,  if  any,  high  contrast  points  and  the  centroid 
algorithm  often  performs  poorly.  Additionally  the  centroid  algorithm  is  sensitive 
to  noise  such  as  that  due  to  the  random  arrival  of  photons  to  the  sensor,  shot  noise. 
Cross  correlation  algorithms  offer  improved  noise  rejection  properties  as  well  as  better 
performance  over  scenes  where  there  are  few  prominent  high  contrast  points,  such  as 
that  found  with  some  extended  sources.  They  do  have  one  serious  drawback  in  that 
the  computational  burden  due  to  the  two  dimensional  cross  correlations  carried  out 
is  often  heavy  and  renders  the  method  slow  and  intensive. 

The  projection  algorithm  as  derived  in  Ref.  [5]  utilizes  two  separate  sensor  arrays 
to  sense  the  tilt  in  each  of  two  dimensions.  Optical  tilt  over  the  incoming  wavefront 
is  fed  to  each  of  the  two  sensors  via  a  beam  splitting  mechanism  and  produces  a  shift 
in  spot  position  on  the  array.  The  projection  based  algorithm  relies  on  the  data  read 
out  from  the  sensors  in  vector  form  and  uses  the  cross-correlation  technique  to  derive 
tilt  in  each  dimension.  This  vector  readout  method  is  not  new  and  has  been  used 
before  by  MIT  Lincoln  Labs  in  their  SWAT  system  [1],  It  does  depart,  however,  in 
the  method  used  to  estimate  optical  tilt,  the  cross  correlation  of  the  vector  read  out. 

2.5  Description  of  Algorithms 

2.5.1  Centroid  algorithm.  The  centroid  algorithm  used  in  this  report  is 
the  traditional  algorithm  whereby  the  center  of  mass  is  calculated  by  dividing  the 
sum  across  both  x  and  y  dimensions  of  the  pixel  values,  weighted  according  to  their 
displacement  from  the  center,  by  the  unweighted  sum  across  both  dimensions.  The 
calculation  is  similar  for  both  the  x  and  y  dimensions.  The  equation  for  the  x  dimen¬ 
sion  is 


Y.xY'vifay) 
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(2.26) 


where  x  is  the  position  in  the  detector  array  from  the  origin  and  i(x,y)  is  the  signal 
at  detector  coordinates^,  y).  As  used  in  simulation  described  in  Sect.  3.2,  this  tilt  is 
in  units  of  detector  pixels.  To  enable  comparison  with  other  algorithms  we  convert  it 
into  rad/meter  in  the  following  manner. 

If  we  take  the  PSF  for  no  tilt  to  be 

2 

(2.27) 


Psf(n)  = 


N 


-j  2irnk 


N 


k=l 


where  A  is  the  amplitude  of  the  kth  pixel  in  the  array,  N  is  the  dimension  of  the  array 
in  pixels  and  n  is  the  pixel  index,  then  the  PSF  for  the  motion  of  one  pixel  is 


Psf(n ) 


N 

E.  ,  ,  .  -j27r(ra-l)fc 

A(k)e  * 


(2.28) 


which  can  be  separated  in  the  exponential  to  give  the  slope  of  the  tilt  change  to 
be  2tt/N  rad/aperture  pixel.  This  is  easily  adapted  to  give  the  conversion  from 
detector  pixels  to  rad/m  as  2nP/(ND)  rad/m/detector  pixel  of  motion,  where  D  is 
the  aperture  diameter  in  meters  and  P  is  in  aperture  pixels. 

The  centroid  algorithm  is  the  optimal  tilt  estimator  for  Gaussian  focal  spots 
with  Poisson  noise  [6].  However  as  can  be  seen  from  Eq.  (2.26),  the  outlying  pixels 
are  weighted  more  than  the  ones  closest  to  the  center.  The  dimmest  pixels  are  often 
found  furthest  from  the  center,  which  means  that  the  least  important  and  noisiest 
pixels  are  given  a  higher  weighting  in  this  algorithm  which  contributes  to  its  lower 
performance  in  noisy  extended  scenes. 


2.5.2  Projection  algorithm.  The  key  advantage  to  using  the  projection- 
based  cross-correlation  algorithm  is  that  the  two-dimensional  image  is  reduced  to  a 
one-dimensional  image,  which  preserves  the  tilt  information  present  in  one  dimension 
of  the  original  image.  So  that  tilt  information  in  both  dimensions  can  be  measured, 
two  sensors  are  required.  The  operation  of  both  is  identical,  and  for  the  purposes  of 
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Figure  2.16:  Projection  transformation  employed  on  a  sensor 

this  report,  only  analysis  for  one  sensor  is  shown.  Referring  to  Figure  2.16,  assuming 
a  W  x  W  lenslet  array  forming  images  onto  a  Z  x  Z- pixel  sensor  array  and  considering 
just  one  image  formed  on  the  sensor,  the  data  are  read  off  in  vector  form  such  that  each 
projection  drk(s)  is  the  summation  of  the  image  across  the  columns  of  the  sub-array 
being  considered  as 


z/w 

<(s)  =  ££>„(  r,s)  (2.29) 

r= 1 

where  (r,  s)  are  sensor  array  coordinates,  Dk(r,s)  is  image  associated  with  the  kth 
observation,  and  Z  and  W  are  as  previously  defined. 

In  Ref.  [6]  it  is  shown  that  the  maximum-likelihood  (ML)  estimator  is  a  good 
slope  estimator  for  wavefronts  in  the  presence  of  noise.  This  combined  with  the 
benefits  of  using  cross-correlation  make  it  beneficial  to  use  a  Bayesian  estimator  for 
the  tilt  parameter.  Here  the  estimator  is  derived  by  forming  the  likelihood  function  for 
the  tilt  conditional  on  the  data  and  then  maximizing  it  with  respect  to  the  tilt  [27]. 
The  likelihood  function  can  be  cast  as  the  conditional  a  posteriori  density  or  the 
conditional  probability  of  the  tilt  parameter,  given  the  measured  projection  and  the 
tilt  parameter  from  the  previous  image  frame.  It  is  assumed  that  an  estimate  of  the 
true  projection  is  known  and  is  deterministic,  and  the  associated  tilt  parameter  (3  in 
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this  case  for  vertical  shift  is  random.  The  estimate  is  given  by 


z/w 

ir(s)wr(s )  =  i(r,  s)w(r,  s )  (2.30) 

r=l 


where  i(r,  s)  is  a  reference  frame,  which  may  be  the  first  frame  captured,  and  wr(s), 
w{r ,  s)  are  window  functions  of  value  1  in  the  area  of  interest  and  zero  otherwise. 

In  order  to  capture  the  discrete  and  non-negative  nature  of  the  photon  counting 
process,  a  Poisson  model  is  used  for  the  measured  projection  dk(s)  corresponding  to 
the  kth  frame  captured.  The  model  having  a  point-wise  mean  of  ir(s  —  (5k)wr{s)  for  a 
given  shift  in  the  vertical  of  (5k-  As  previously  stated,  the  likelihood  function  requires 
knowledge  of  the  distribution  of  the  tilt  parameter  from  frame  to  frame.  This  can  be 
illustrated  using  Bayes’  rule,  in  that 


fPk\drk,pk-i(P\dib') 


fd-l  (d) 


(2.31) 


where  f/3k\drk,/3k_1(b\d,  b')  is  the  probability  of  the  tilt  for  frame  k  conditioned  on  the 
measured  projection  drk  and  the  tilt  from  the  previous  frame  Pk-i  ,  f f3k\Pk-\(P\b')  is  the 
probability  of  the  random  tilt  being  equal  to  (5k  =  b,  conditioned  on  the  tilt  parameter 
of  the  previous  frame  being  (5k-\  =  b' ,  fdrk\gk(d\b)  is  the  probability  that  the  projection 
vector  random  process  is  equal  to  a  specific  realization  of  that  process  conditioned  on 
(5k  =  b  and  fdrk(d )  is  the  unconditional  probability  that  the  projection  vector  random 
process  is  equal  to  a  specific  realization. 

To  estimate  the  tilt  parameter,  Eq.  (2.31)  is  maximized  with  respect  to  its 
logarithm  L{b)  giving  the  maximum  likelihood  function 


L(b)  =  \n[fdr]Pk(d,b)]  +  ln[f pk\pk  l(b\b')\  -  In [fdl(d)],  (2.32) 
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which  can  be  expressed  as  a  function  of  terms  that  only  depends  on  the  shift  b,  by 
dropping  the  third  term  which  is  independent  of  b  to  give 


L{b)  =  In  [f<rk\/}k(d>b)\  +  Hfpk\Pk-i(b\b%  (2-33) 


In  order  to  maximize  this  expression  in  terms  of  the  tilt  parameter  b,  we  must  have 
some  knowledge  of  the  probability  of  the  current  tilt  parameter  conditioned  on  the 
tilt  from  the  previous  frame.  If  this  knowledge  is  available  it  should  be  used,  however 
in  many  cases  it  is  not  and  it  is  common  practice  [6]  to  choose  a  uniform  density 
which  does  not  vary  with  b  and  the  second  term  of  Eq.  (2.33)  may  be  dropped.  In 
this  research  a  uniform  window  centered  on  the  previous  tilt  estimate  is  used. 

The  projection  vector  dk  represents  an  ensemble  of  independent  Poisson  random 
variables  associated  with  individual  pixel  measurements.  Given  the  initial  assumption 
of  statistical  independence  between  measurements,  the  pdf  of  a  collection  of  samples 
of  drk  given  (ik  =  b  can  be  expressed  as  a  product  of  the  marginal  densities  over  all 
pixels  in  dk,  such  that 


z/w 


/»(<*.(>)  =  P«  =  <i|  0„  =  b)=  JI 


w  (sh  is 


—  b)d^e  ^ 


(s—b)wr(s) 


s= 1 


d(s)\ 


(2.34) 


giving  the  log-likelihood  function  to  be 


z/w 

L(b)  =  ^  d(s)  ln[icr(s)ir(s  —  b )]  —  wr(s)ir(s  —  b )  (2.35) 

S=1 


In  this  research  the  widowing  function  wr(s )  is  chosen  to  be  smaller  than  the 
size  of  the  projection  vector  by  a  number  of  pixels  equal  to  2 5max,  where  bmax  is  the 
expected  maximum  absolute  value  of  the  tilt  parameter.  Then,  because  the  window 
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function  is  defined  as 


wr(s ) 


0  for  1  <  s  <  6max, 

<  0  for  Z/W  -  bmax  <s<  Z/W , 
1  elsewhere, 


the  natural  logarithm  is  equal  to  negative  infinity  for  values  of  s  which  make  the 
window  function  equal  to  0.  To  avoid  this,  the  limits  of  integration  are  chosen  to 
include  only  those  points  where  the  window  function  is  non-zero,  making  the  log- 
likelihood  function 


Z  /W  —  femax 

m=  E  d/s)  ln[ir(s  —  b)]  —  ir(s  —  b).  (2.36) 

5— 6max 


Equation  (2.36)  can  be  maximized  with  respect  to  tilt  parameter  b  using  an 
iterative  approach  to  computing  the  value  of  the  function  locally  around  the  current 
estimate  for  b  and  updating  the  estimate  in  the  direction  of  increasing  L  in  steps  of 
A b.  The  value  of  A b  becomes  the  resolution  of  the  tilt  estimation  algorithm  in  units 
of  array  pixels.  A  linear  interpolator  is  chosen  to  produce  sub-pixel  resolution  for  tilt 
estimates.  The  linear  interpolator  has  the  form 

ir(s  -  b)  =  [(1  -  bf)ir(s  -  bi )  +  (■ bf)ir(s  -  (6*  -  1))],  (2.37) 

where  bi  and  bf  are  the  integer  and  fractional  parts  of  b  respectively. 

In  general,  the  wavefront  sensing  and  tilt  estimation  problem  has  been  addressed 
well  in  the  case  of  point  sources  and  to  some  extent  extended  beacons  [20,21],  The 
projection-based  cross  correlation  algorithm  has  been  suggested  as  being  effective  in 
viewing  stellar  and  laser  beacons  but  little  has  been  done  to  evaluate  its  performance 
in  the  specific  case  of  speckle  noise.  Research  has  been  done  on  the  effects  of  speckle 
on  various  aspects  of  motion  estimation  and  imaging  Ref.  [2,7,25,28].  This  research 
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aims  to  examine  that  particular  case  where  the  projection  algorithm  is  applied  to  the 
speckle  noise  case. 
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III.  Simulation  Structure  and  Models 


This  chapter  conveys  the  approach  and  methodology  used  to  address  the  research 
goals.  First,  the  simulation  structure  is  presented  and  discussed.  Then,  each 
conceptual  model  used  in  the  simulation  is  explained.  Lastly,  the  metrics  developed 
to  enable  meaningful  conclusions  to  be  drawn  from  the  results  are  presented  and 
discussed. 

3. 1  Simulation  Structure 

This  section  describes  how  the  simulation  is  built  and  provides  a  visual  repre¬ 
sentation  of  the  physical  simulation  structure  to  assist  the  reader  in  conceptualizing 
the  simulation  goals.  The  simulation  goals  are  to  compare  and  contrast  the  perfor¬ 
mance  of  projection-based  cross-correlation  and  traditional  centroiding  tilt  estimation 
algorithms  in  the  presence  of  speckle  noise.  Figures  3.1  and  3.2  show  the  conceptual 
structure  of  the  simulation  for  each  of  the  centroiding  and  projection-based  cases. 
Simply,  the  simulation  involves  a  source  propagated  to  an  optical  element  which  fo¬ 
cusses  the  light  onto  a  wavefront  sensor.  In  both  cases,  the  source  is  situated  on  the 
optical  axis,  a  distance  Z  away  from  the  focussing  element,  a  thin  lens.  It  is  assumed 
that  Z  3>  2D2/X  where  D  is  the  diameter  of  the  pupil  in  meters  and  A  is  the  wave¬ 
length  of  light  in  question.  This  geometry  was  chosen  to  simplify  the  propagation 
problem  to  one  of  simply  performing  a  Fourier  transform  to  propagate  the  source  to 
the  lens.  It  was  not  a  primary  goal  of  this  research  to  compare  the  results  with  and 
without  turbulence.  However,  real-world  scenarios  would  always  have  some  degree  of 
atmospheric  turbulence  present,  and  for  this  reason  it  is  given  consideration  here.  The 
turbulence  in  this  simulation  is  applied  at  the  lens  surface.  The  turbulence  applied 
is  assumed  to  simulate  path  turbulence  between  the  source  and  the  lens.  For  simu¬ 
lation  purposes  the  wavefront  sensor  is  located  at  the  focal  plane  of  the  lens.  Again, 
this  allowed  a  Fourier  transform  to  be  used  to  propagate  from  the  lens  to  the  sensor. 
In  the  case  of  the  projection-based  algorithm  shown  in  Figure  3.2,  a  beam-splitter 
is  situated  between  the  lens  and  each  wavefront  sensor.  The  reasoning  for  splitting 
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Figure  3.1:  Conceptual  Diagram  of  Simulation  for  the  Centroiding  case 


Table  3.1:  Simulation  Parameter  Space 


Parameter 

Symbol 

Units 

Value 

Propagation  grid 

N 

pixels 

128 

Aperture  Dimension 

P 

pixels 

64 

Aperture  Diameter 

D 

m 

0.07 

Fried  Seeing  Parameter 

r0 

m 

0.07 

Intensity  DOF 

M 

unit  less 

1 

Window  function 

w 

pixels 

24 

the  beam  into  separate  X  and  Y  plane  components  is  described  in  Section  2.5.2.  In 
this  simulation  only  one  sensor  is  simulated  as  the  performance  of  each  is  identical. 
Table  3.1  lists  the  parameters  used  throughout  the  simulation. 

3.1.1  Simulation  Description.  This  section  describes  the  process  of  the 
simulation.  Both  the  projection-based  and  centroiding  algorithms  are  simulated  in 
each  of  the  following  four  scenarios: 


1.  Extended  source  without  turbulence, 

2.  Speckle  source  without  turbulence, 
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Figure  3.2:  Conceptual  Diagram  of  Simulation  for  the  Projection-based  case. 

3.  Extended  source  with  turbulence,  and 

4.  Speckle  source  with  turbulence. 

The  analysis  computes  the  rms  tilt  error  for  mean  light  levels  of  100  to  1000  photons 
received,  in  steps  of  100  photons.  At  each  light  level,  1000  realizations  of  the  propa¬ 
gation  are  used  to  compute  the  rms  tilt  error.  As  the  tilt  variance  can  be  very  large, 
the  rms  error  is  computed  to  better  model  the  tilt  measured  by  the  sensor  due  to 
pixel  integration  time.  A  separate  phase  screen,  simulating  atmospheric  turbulence, 
is  used  for  each  realization.  This  increases  the  randomness  of  the  simulation.  Poisson- 
distributed  random  noise  is  added  to  each  realization  to  simulate  the  random  arrival 
nature  of  the  photons,  known  as  shot  noise.  A  non-noisy  image  frame  is  used  as  a 
reference  at  each  light  level.  For  the  projection-based  model,  the  reference  frame  is 
computed  from  the  average  of  1000  realizations  of  the  source  to  image  plane.  This 
corresponds  the  the  idea  that  the  sensor  may  be  staring  at  the  target  for  a  period  of 
time  before  tracking  commences.  The  original  source  image  is  used  for  the  centroid- 
ing  algorithm,  since  the  source  is  on  the  optical  axis.  Then  the  tilt  measured  by  the 
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centroid  on  the  source  image  is  the  zero  tilt  reference  for  the  centroiding  algorithm. 
Listing  III.  1  shows  the  Matlab®  code  used  to  compute  the  reference  image. 

Listing  III.  1:  Reference  Image  Code  Section.  (chapter3/ReferenceImageCode.m) 

1  otfl  =  otf  .array  (  :  ,  :  ,  iloopl )  ;  °/n  choose  otf  from  tilt  removed  phase... 

screen 

if  iloopl  ==  1 

for  iloop3  =  1: images 

psfl  =  ( abs ( if f 1 2 ( ( f f 1 2 ( ext  ended. sour ce ) )  . *  ... 

otfl)))  ; 

psf  l_correct_kbar  =  psfl  *  Kbar  (  iloop2  )  ;  °/»  set... 

average  rxd  photon  count  Kbar 

6  temp=ones (128,128) ,/(ones(128,128)+... 

psf l_correct_kbar /2)  ; 

speck_image_correct_Kbar  =  icdf (  ’  nbin  ’  , rand . .  . 

(128,128)  ,1,  temp)  ;  set  average  rxd  photon... 
count  Kbar 

speckle.ref .image  =  speckle.ref .image  +  (... 

speck_image_correct_Kbar) / images ; 
pro j _ref .image  =  pr oj _ref _ image  +  (... 
psf l_correct_ kbar)/ images ; 

11  end  1  for  iloop3 

end  7.  if  iloopl 


In  order  to  propagate  the  source  through  the  system,  firstly,  the  Matlab® 
function  Make-otf.m,  described  in  Section  3.2.1  is  used  to  compute  the  Optical  Trans¬ 
fer  Function  (OTF)  of  the  propagation.  Then  the  inverse  Fourier  transform  of  the 
product  of  the  source  and  the  OTF  is  computed  to  give  the  PSF  which  is  normalized 
for  one  photon.  At  this  point  the  PSF  is  corrected  for  the  desired  mean  light  level 
(100  to  1000  photons).  The  speckle  image  is  now  created  as  described  in  Section  3.2.3. 
Shot  noise  is  then  applied  using  the  Matlab®  function  poissrnd.  The  wavefront  tilt 
estimates  using  both  centroid  and  projection-based  methods  are  then  calculated  using 
the  functions  centroid. m  and  projection_methodX.m,  which  are  described  in  Sections 
3.2.5  and  3.2.6,  respectively.  The  tilt  error,  defined  here  as  the  difference  between  the 
known  or  reference  tilt  and  the  computed  tilt  of  the  new  image,  is  then  computed  for 
each  scenario.  The  projection-based  algorithm  internally  computes  the  tilt  in  rad/m. 
However,  the  centroiding  algorithm  does  not.  Therefore,  the  measured  tilt  in  the 
centroiding  case  must  be  converted  into  units  of  rad/m  before  the  error  from  both 
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methods  is  stored  and  used  to  compute  the  rms  tilt  error  as  described  in  Section  3.3. 
The  rms  tilt  error  for  each  light  level  is  then  stored  for  display  in  plot  form. 

3.2  Simulation  Models 

3.2.1  Optical  Transfer  Function  Model.  The  system  OTF  is  computed  by 
the  function  Make.otf2.rn  which  is  shown  in  Listing  A. 2.  Make_otf2.m  works  by 
creating  an  aperture  function  and  filling  it  with  the  phase  at  the  aperture  (provided 
by  the  user).  In  this  simulation,  the  phase  is  the  turbulence- induced  wavefront  error 
produced  by  the  phase  screen  generation  code  as  described  in  Section  3.2.4  or  in 
the  case  where  turbulence  is  not  considered,  a  zero  phase  screen.  The  function  then 
performs  an  autocorrelation  of  the  pupil  using  the  fft,2  function  and  then  normalizes 
the  result  to  yield  the  OTF. 

3.2.2  Extended  Source  Model.  An  often-used  definition  of  an  extended 
source  is  one  where  the  source  is  considered  extended  when  it  can  be  resolved  in  the 
image  plane.  A  considerably  more  in-depth  description  is  provided  in  Roggemann  and 
Welsh  [21],  and  further  references  to  extended  sources  are  made  in  Refs.  [12,17,19, 
20,22],  Although,  none  of  those  cited  works  gives  clear  guidance  on  when  a  source  is 
extended  and  when  it  is  not.  For  the  purposes  of  this  research,  all  that  is  required  is 
for  the  source  to  be  sufficiently  large  to  adequately  display  the  effects  of  speckle  noise 
and  to  be  resolved  in  the  image.  For  this  reason,  the  extended  source  was  chose  to  be 
a  4  x  4  pixel  source  in  a  128  x  128  propagation  grid  which,  when  imaged  through  the 
system,  was  resolved  and  provided  an  acceptable  level  of  speckle  distortion.  Figures 

3.3  and  3.4  show  the  extended  source  model  and  the  image  produced  by  the  simulation. 
It  clearly  meets  the  objective  of  being  resolved. 

3.2.3  Speckle  Noise  Model.  Various  statistical  models  for  use  in  simulating 
speckle  noise  are  described  in  Section  2.3.3.  The  negative  binomial  model  was  chosen 
for  this  simulation  because  of  its  relative  ease  of  implementation  in  Matlab®  .  No 
direct  function  exists  for  generating  negative  binomial  distributed  random  variables  in 
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Figure  3.3:  Extended  source  Model 
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Figure  3.4:  Image  of  extended  source 
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Matlab®  .  However,  there  is  an  Inverse  Cumulative  Distribution  Function  ( ICDF ), 
which,  when  given  an  appropriate  probability  distribution  and  other  parameters,  gen¬ 
erates  a  random  variable  of  the  desired  distribution  and  mean.  The  ICDF  returns  an 
array  of  values  of  the  inverse  cumulative  distribution  for  a  specified  probability  distri¬ 
bution,  given  an  input  data  set  and  mean.  In  this  simulation,  the  desired  distribution 
is  Negative  Binomial,  the  data  set  is  a  random  realization  of  128  x  128  pixel  image 
with  pixel  intensity  from  0  to  1  and  the  mean  is  a  128  x  128  array  in  which  each  value 
is  set  according  to 


r>  M 
M  +  K ’ 


(3.1) 


where  B  is  the  mean  and  A4  and  K  are  as  previously  defined. 


This  method  of  generating  random  variables  has  been  validated  through  numer¬ 
ical  analysis  [13].  A  128  x  128  pixel  image  (later  windowed  to  24  x  24  as  described  in 
3.2.6)  was  computed  by  propagating  the  source  using  a  series  of  Fourier  transforms 
and  applying  an  OTF  generated  from  the  OTF  model  from  Section  3.2.1  to  compute 
the  PSF.  The  PSF  is  then  corrected  for  the  desired  mean  light  level.  A  temporary 
array  is  created  with  the  correct  mean  photon  count  and  supplied  to  the  icdf  func¬ 
tion  to  produce  the  speckle  image  with  correct  statistical  distribution  and  mean  light 
level.  Listing  III. 2,  taken  from  the  parent  hie  Thesis  simulation  _code  Ml  1,  is  used  to 
implement  this  model.  Figure  3.5  shows  an  example  of  the  image  of  a  speckle  source 
created  with  the  previously  discussed  method. 

Listing  III. 2:  Speckle  Model  Code  section.  (chapter3/Specklemodel.m) 

psfl  =  ( abs ( if f 1 2 ( ( f f 1 2 ( ext  ended. sour ce ) )  . *  otfl))); 

3  psf  l_correct_kbar  =  psfl  *  Kbar  ( iloop2 )  ;  °/0  set  average  rxd... 

photon  count  Kbar 

t emp  =  one  s (128,128)  ./(ones(128,128)+psfl_correct_kbar/2)  ; 

speck_image_correct_Kbar  =  icdf(,nbin,,rand(128,128),l,... 
temp);  °/n  set  average  rxd  photon  count  Kbar 


3.2.4  Atmospheric  Turbulence  Model.  Atmospheric  turbulence  results  in 
variations  in  the  refractive  index  along  optical  path.  This  variation  is  a  random 
process,  and  so  a  model  of  this  turbulence  should  be  a  statistical  average  of  this 
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Figure  3.5:  Image  of  speckle  source 

random  variation.  Phase  screen  creation  involves  generating  single  realizations  of  this 
random  process.  This  can  be  achieved  by  transforming  computer-generated  random 
numbers  into  two  dimensional  arrays  of  phase  values  that  have  the  same  statistical 
distribution  as  the  random  turbulence-induced  phase  variations.  In  most  cases  the 
phase  is  written  as  a  sum  of  basis  functions.  Two  commonly  used  basis  sets  are 
Zernike  polynomials  and  Fourier  series.  The  Fourier  series  is  the  basis  set  used  in  this 
simulation.  A  detailed  description  of  this  method  is  found  in  Ref.  [11],  Briefly,  this 
method  involves  writing  the  optical  phase  4>(x,  y )  as  a  Fourier  series: 

OO  OO 

4>{x,y)  =  X  X  mexp[i27r(fXrix  +  fVmy)\,  (3.2) 

n= — oo  m=— oo 

where  fXn  and  fVm  are  the  x  and  y  spatial  frequencies  and  cn,m  are  the  Fourier  coeffi¬ 
cients.  Treating  the  phase  as  a  two  dimensional  signal  and  making  use  of  Parseval’s 
theorem  the  Fourier  series  coefficients  become 

(|c„,m|2>  =  (3.3) 

where  &(fXn,  fym)  is  the  power  spectral  density  of  the  turbulence-induced  phase  delay 
and  A fXn  and  A fym  are  the  corresponding  sample  spacings  of  the  spatial  frequen¬ 
cies.  The  expectation  has  been  taken  because  the  phase  is  a  random  process.  Given 
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the  phase  variations  are  independent  from  one  another  in  a  statistical  sense,  the 
central-limit  theorem  can  be  used  to  determine  that  the  coefficients  are  Gaussian  in 
distribution,  which  obey  circular  Gaussian  statistics  of  zero  mean  and  variance  given 
in  Eq.  (3.3).  Knowing  this,  random  Gaussian  numbers  generated  via  Matlab®  can 
be  multiplied  by  the  square  root  of  the  variance  in  Eq.  (3.3).  In  this  simulation  the 
Fourier-series  coefficients  are  generated  on  a  uniformly  sampled  grid  and  the  FFT 
method  is  used  to  synthesize  phase  screens  with  a  good  degree  of  computational  ef¬ 
ficiency.  Listing  A. 7  gives  the  Matlab®  code  used  to  implement  this  method  and 
generate  the  phase  screens.  However  as  noted  in  Ref.  [11]  this  method  can  result 
in  poor  simulation  of  low-order  modes  such  as  tilt.  In  this  research  it  is  desired  to 
only  examine  the  tilt  induced  by  speckle  noise,  and  not  introduce  anymore  tilt  than 
that  of  the  speckle.  Therefore  this  is  not  an  issue.  Furthermore,  any  residual  tilt  is 
removed  from  the  phase  screen  by  generating  an  array  with  the  projected  tilt  across 
it  and  subtracting  this  from  the  generated  phase  screen.  Listing  III.3  from  Listing 
A. 2  shows  Matlab®  code  used  to  implement  the  tilt  removal. 

Listing  III. 3:  Tilt  Removal  Code  section.  (chapter3/TiltRemovalCode.m) 

7.7.  tilt  removal  section 
tx  =  2*pi*(-64:63) / 128 ; 
tx_matrix  =  one s ( 1 28 , 1 ) * tx ; 

wvs_x  =  sum ( sum ( tx_matr ix . *  phase .* aperture ))  ; 

5  wvs_x  =  wvs_x/ sum ( sum ( aperture .* tx.matrix . " 2 )) ; 
new_screen  =  phase  -  wvs_x*tx_matrix ; 

Figure  3.6  shows  an  example  phase  screen  generated  by  this  model. 

3.2.5  Centroid  Method.  The  centroiding  algorithm  described  in  Section 
2.5.1  is  implemented  in  MATLAB®  code.  The  centroiding  equation,  Eq.  (2.26)  is 
almost  coded  line  for  line  in  Matlab®  .  The  supplied  image  is  windowed  in  the 
same  dimensions  as  the  projection-based  method  uses  to  enable  valid  comparison  of 
results.  Listing  A. 3  shows  the  Matlab®  code  used.  The  code  outputs  the  tilt  in 
units  of  pixels  so  this  must  be  converted  to  rad/m,  as  described  in  Section  3.3,  for 
comparison  with  the  projection  result.  As  in  the  projection-based  method,  only  one 
dimension  is  simulated  here. 
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Figure  3.6:  Sample  phase  screen  created  by  the  FFT  method  in  this  simulation. 

3.2.6  Projection-based  Method.  There  are  two  functions  in  this  simulation 
that  implement  the  projection  algorithm  as  described  in  Section  2.5.2.  They  are 
functionally  identical,  differing  only  in  the  manner  in  which  the  linear  interpolation 
step  is  implemented,  the  reasons  for  this  are  explained  later  in  this  section.  The 
algorithm  is  implemented  in  two  steps.  Step  one  is  the  formation  of  the  vector  drk(s) 
as  in  Eq.  (2.29),  by  summing  across  the  sensor  in  the  dimension  desired,  x  or  y, 
which  would  correspond  to  the  sensor  being  read  out  in  vector  mode.  Lines  9  to 
11  of  Listings  A. 4  and  A. 5  show  how  this  is  implemented.  This  process  is  carried 
out  on  both  the  image  frame  of  interest  and  the  reference  image  frame.  Step  two  is 
the  vector  cross-correlation,  which  involves  the  formation  of  the  likelihood  function 
and  the  linear  interpolation  as  described  in  Section  2.5.2  and  Eq.  (2.33),  (2.34)  and 
(2.35).  The  windowing  of  the  log-likelihood  function  as  described  in  Section  2.5.2  and 
shown  in  Eq.  (2.36)  is  implemented  by  line  8  and  the  loop  beginning  at  line  14  of 
listings  A. 4  and  A. 5.  Note  the  end  points  and  step  size  of  the  loop  are  identical  to 
those  of  the  vector  initialized  in  line  8.  These  lines  of  code  set  the  window  size  to 
±5  pixels  from  image  center,  which  corresponds  to  the  optical  axis.  The  window  size 
is  modeled  on  that  chosen  in  the  proposed  implementation  of  this  algorithm  in  the 
Dunn  Solar  Telescope  in  Ref.  [5].  In  that  case,  an  E2V  CCD60  array  with  pixels  of 
24 /xm  pitch,  and  a  sub  array  size  of  24  x  24  pixels  was  used,  yielding  a  held  of  view 
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of  1/3  of  an  arc  second.  The  difference  in  the  two  functions  came  from  experience 
gained  during  the  simulation  evaluation.  The  difference  lies  in  the  method  used  to 
derive  sub-pixel  accuracy  in  the  linear  interpolation.  Both  methods  achieve  the  same 
accuracy  using  different  approaches.  ProjectionMethod2  separates  the  interpolation 
into  an  integer  and  fractional  component  and  interpolates  from  there.  This  achieved 
good  results  with  the  speckle  case,  and  poorer  results  with  the  non-speckle  case. 
Through  research  it  was  found  that  inherent  sub-pixel  accuracy  of  the  makeshiftjvec 
function  achieved  better  results  with  the  non-speckle  case  and  thus  two  versions  are 
used. 


3.3  Simulation  Metrics 

In  order  to  make  proper  comparison  of  results  obtained,  the  rms  tilt  error  due 
to  noise  in  units  of  radians  per  meter  was  calculated  for  all  simulations.  As  previously 
mentioned  the  projection-based  method  inherently  computes  the  tilt  error  in  units  of 
rad/m.  By  examining  the  centroiding  equation,  Eq.  (2.26),  it  can  be  seen  that  the 
tilt  computed  by  the  centroiding  method  is  in  units  of  detector  pixels.  To  enable 
comparison  with  other  algorithms  we  convert  it  into  rad/m  in  the  following  manner. 
If  we  take  the  point  spread  function  (PSF)  for  no  tilt  to  be  as  shown  in  (2.27)  then 
the  PSF  for  motion  of  one  pixel  is 


PS F(n  -  1) 


E  -jSn(n-l)k 

A(k)e  N 


k= 1 


2 


5 


(3.4) 


which  can  be  separated  in  the  exponential  to  give  the  slope  of  the  tilt  change  to 
be  2n/N  rad/aperture  pixel.  This  is  easily  adapted  to  give  the  conversion  from 
detector  pixels  to  rad/m  as  2nP/(ND)  rad/m/Detector  pixel  of  motion,  where  D  is 
the  aperture  diameter  in  meters  and  P  is  in  aperture  pixels. 
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The  rms  tilt  error  an  is  computed  by 


(3.5) 


where  Nt  is  the  number  of  realizations  made  at  each  mean  light  level  and  A  is  the 
computed  tilt  error  for  a  given  realization. 
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IV.  Simulation  Results 


This  chapter  presents  the  results  of  the  simulation  described  in  Chapter  III  and 
presents  some  discussion  in  the  context  of  the  original  research  goals. 

4-1  Simulation  using  Extended  Source  No  Speckle 

The  research  done  by  Cain  in  Ref.  [5]  and  Cain,  Hyatt,  and  Armstrong  in 
Ref.  [4]  established  a  baseline  premise  for  the  results  of  this  research,  in  that  the 
projection-based  cross-correlation  algorithm  exhibited  improved  performance  when 
viewing  stellar  beacons  as  well  as  extended  objects  and  scenes.  The  first  goal  of  this 
research  was  to  establish  the  performance  of  the  projection-based  algorithm  when 
considering  extended  objects.  Figure  4.3  shows  the  computed  tilt  error  for  both 
the  projection-based  and  traditional  centroiding  algorithms  in  plot  (a)  the  case  with 
no  turbulence  considered  and  plot  (b)  with  turbulence  considered.  In  Figure  4.3(a) 
it  can  be  seen  that,  as  expected,  the  projection  algorithm  performs  better  than  the 
centroiding  algorithm  at  all  mean  light  levels  simulated.  Although  the  difference  is  not 
significant,  it  is  measurable.  In  Figure  4.3(b)  it  can  be  seen  that,  again  the  projection 
method  exhibits  better  performance,  in  this  case  by  an  increased  margin.  Whilst  the 
computed  error  has  increased  for  both  methods,  the  increase  in  error  produced  by  the 
centroid  method  was  much  larger  and  tends  to  stabilize  as  the  light  level  increases 
while  the  project  method  continues  to  improve. 

4-2  Simulation  using  Extended  Source  with  Speckle 

Figure  4.2  shows  the  computed  tilt  error  for  both  the  projection-based  and 
traditional  centroiding  algorithms  in  the  presence  of  speckle  noise  for  plot  (a)  the 
case  with  no  turbulence  considered  and  plot  (b)  turbulence  considered.  Again,  the 
projection  method  performs  better  than  the  centroid,  although  with  a  much  less 
observable  difference.  This  indicates  that  the  speckle  noise  is  the  dominant  effect, 
and  that  the  turbulence  may  even  exacerbate  the  effect  of  speckle. 
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Figure  4.1:  (a)  Tilt  error  for  extended  source  with  no  turbulence, 

(b)  Tilt  error  for  extended  source  with  turbulence. 
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Figure  4.2:  (a)  Tilt  error  for  speckle  with  no  turbulence, 

(b)  Tilt  error  for  speckle  with  turbulence. 
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Figure  4.3:  (a)  No  turbulence  considered  tilt  error. 

(b)  Turbulence  considered  tilt  error. 

4-3  Effect  of  Atmospheric  Turbulence  and  Speckle 

When  considering  all  results,  it  can  be  seen  that  once  speckle  is  introduced  the 
turbulence  effect  is  minimal  by  comparison.  The  projection-based  algorithm  consis¬ 
tently  performs  better  than  the  traditional  centroiding  algorithm  in  all  cases.  It  is 
interesting  to  note  however,  that  the  increasing  light  levels  did  not  appear  to  favor 
either  method  when  speckle  was  considered. 
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V.  Conclusions  and  Further  Work 


This  chapter  presents  a  summary  of  the  research  and  key  results  from  Chapter  IV. 

Further,  it  addresses  the  main  challenges  faced  during  this  research  and  how 
well  the  initial  goals  were  satisfied.  Recommendations  of  further  work  that  could 
begin  with  and  extend  on  these  results  is  also  discussed. 

5. 1  Summary 

The  objective  of  this  research  was  to  examine  the  performance  of  a  projection- 
based  cross  correlation  algorithm  as  developed  by  Cain  in  Ref.  [5]  in  the  presence  of 
speckle  noise  via  simulation.  Three  areas  of  research  made  up  the  bulk  of  this  work. 
They  were: 

1.  Extended  Beacon  modeling, 

2.  Speckle  noise  modeling,  and 

3.  Implementation  of  MLE  and  the  projection-based  algorithm. 

It  was  necessary  to  have  all  three  areas  functioning  correctly  to  produce  results 
from  which  valid  conclusions  could  be  drawn.  As  discussed  in  Chapter  III,  consid¬ 
eration  was  also  given  to  effect  of  atmospheric  turbulence  on  the  performance  of  the 
algorithm. 

The  research  began  with  an  examination  of  the  performance  of  the  projection 
algorithm  with  extended  beacons.  This  required  the  development  of  a  suitable  model 
for  use  in  simulation.  Much  work  [22,28,29],  has  been  done  on  the  theory  of  extended 
versus  point  sources,  but  little  has  been  done  on  mathematically  modeling  them.  This 
presented  a  challenge  in  this  research.  It  was  eventually  decided  that  any  source  that 
could  be  resolved  in  the  final  image  would  be  sufficient  to  model  an  extended  source 
for  the  purposes  of  this  research.  The  performance  of  the  projection  algorithm  was 
examined  in  simulation  against  the  performance  of  the  traditional  centroiding  algo¬ 
rithm  and  was  shown  to  offer  measurable  improved  performance  over  the  centroiding 


51 


method  in  both  the  case  where  turbulence  is  not  considered  and  when  it  was  con¬ 
sidered.  In  particular  the  performance  improvement  margin  was  significant  in  the 
case  where  turbulence  was  considered.  This  could  be  attributed  to  the  ability  of  the 
projection  based  algorithm  to  reject  the  low  spatial  frequency  changes  in  the  image 
introduced  by  turbulence.  Also  noteworthy  in  the  case  where  turbulence  is  consid¬ 
ered,  is  the  reducing  error  trend  as  the  mean  light  level  increases  for  the  projection 
method.  This  contrasts  with  result  for  the  centroiding  method,  which  is  trending  to 
level  out. 

Considerable  work  had  been  done  previously  on  both  characterizing  the  statis¬ 
tical  nature  of  the  speckle  phenomena  [2,8-10],  and  on  the  use  of  MLE  for  motion 
estimation  [3,7,25].  Little  work  had  been  done  on  combining  these  experimentally 
in  simulation  and  examining  the  results.  Several  statistical  models  for  speckle  phe¬ 
nomena  have  been  put  forward  and  in  some  cases  [9, 10]  validated  by  mathematical 
evaluation.  These  models  all  had  limiting  cases  and  applications  and  remained  unclear 
as  to  whether  there  was  an  overall  model  that,  while  not  characterizing  all  aspects 
of  speckle  noise  completely,  would  be  acceptable  for  widespread  use.  Two  models 
in  particular  were  evaluated  for  their  ability  to  model  speckle  in  the  context  of  this 
research,  that  is,  the  case  of  tracking.  The  negative  exponential  distribution  for  pho¬ 
ton  count,  as  discussed  in  Section  2.3  with  later  addition  of  Poisson  distributed  shot 
noise,  was  initially  explored.  However  this  method  produced  images,  Figure  5.1,  after 
propagation  that  were  no  longer  acceptable  as  speckle  images  because  they  did  not 
display  the  known  properties  of  speckle  noise.  The  negative  binomial  distribution 
as  detailed  in  Section  2.3  and  Ref.  [9],  which  includes  both  speckle  effect  and  shot 
noise  was  found  to  generate  images  with  the  correct  speckle  features  as  shown  in 
Figure  3.5.  Using  this  negative  binomial  distribution  reduced  the  computational  load 
in  computing  the  final  image  to  present  to  the  projection  algorithm.  It  was  shown 
that  the  projection  algorithm  again  consistently  performed  better  than  the  traditional 
centroiding  algorithm  in  both  the  case  of  with  turbulence  and  without  turbulence.  It 
is  noteworthy  again  that  the  addition  of  turbulence  had  little  effect  on  the  overall 
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Figure  5.1:  Initial  Speckle  Model  Image 

performance  of  either  algorithm  when  speckle  was  present.  Turbulence  did  not  seem 
to  affect  the  average  performance  of  either  algorithm  measurably,  however  it  did  pro¬ 
duce  a  pronounced  increase  in  the  variation  in  computed  error.  Figure  5.2  shows  the 
results  obtained  with  the  overall  mean  error  for  each  case.  This  would  indicate  that 
speckle  is  the  dominant  effect,  which  is  not  surprising  as  the  statistical  pdf  used  for 
the  MLE  was  Poisson  and  does  not  account  for  the  Bose-Einstein  nature  of  speckle 
noise.  The  reason  for  the  increase  in  error  at  a  mean  light  level  of  >  800  photons  is 
unclear  and  may  in  fact  be  a  statistical  anomaly  which  would  be  resolved  by  more 
realizations  at  each  light  level  with  the  added  burden  of  increase  computation  time 
or  may  be  averaged  out  with  an  increase  in  the  maximum  mean  light  level  simulated 
to  greater  than  1000  photons.  Additional  research  may  be  carried  out  to  determine 
if  there  are  limits  (high  or  low  mean  photon  counts),  outside  of  which  the  existing 
algorithm  no  longer  provides  a  good  approximation.  Overall  the  projection-based 
cross-correlation  algorithm  has  been  shown  to  have  increased  performance  over  the 
traditional  centroiding  algorithm  in  the  presence  of  speckle  noise. 

5.2  Key  Results 

This  section  summarizes  the  key  conclusions  from  the  simulation  results  in 
Chapter  IV. 
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Figure  5.2:  (a)  Computed  mean  tilt  error  no  turbulence, 

(b)  Computed  mean  tilt  error  with  turbulence. 


•  The  projection-based  cross-correlation  algorithm  showed  improved  performance 
in  tilt  estimation  over  the  traditional  centroiding  method  when  considering  ex¬ 
tended  beacons  even  when  atmospheric  turbulence  is  considered  and  the  correct 
(Bose-Einstein)  pdf  was  not  used. 

•  The  projection-based  cross-correlation  algorithm  showed  improved  performance 
in  tilt  estimation  over  the  traditional  centroiding  method  in  the  presence  of 
speckle  noise.  Different  models  for  speckle  noise,  while  statistically  correct, 
may  not  accurately  model  the  phenomena  in  the  case  being  considered. 

•  While  the  projection  algorithm  did  outperform  the  centroid  algorithm  in  all 
cases,  the  difference  in  performance  was  not  always  significant.  With  this  in 
mind,  and  taking  into  account  the  computational  burden  of  the  projection-based 
algorithm,  in  some  cases  depending  on  the  conditions  the  centroiding  algorithm 
may  be  all  that  is  required. 
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5.3  Recommendations 


This  section  provides  recommendations  for  future  work  based  on  this  research. 
Develop  Bose-Einstein  MLE  Algorithm.  Initially  a  goal  of  this  research,  the  MLE 
algorithm  could  be  developed  using  the  Bose-Einstein  statistical  distribution  for  pho¬ 
ton  count  in  speckle  images.  While  this  research  has  shown  that  the  projection-based 
method  has  improved  performance  over  the  traditional  centroiding  method  using  the 
Poisson  distribution,  it  is  not  the  true  distribution  and  the  algorithm  may  show  greater 
improvement  if  the  correct  statistical  model  where  used  in  the  MLE. 

Continuation  of  Speckle  Model.  More  work  could  be  done  characterizing  the 
speckle  noise  model  so  it  could  be  propagated  in  the  same  manner  as  the  extended 
source.  Then  the  algorithms  could  be  further  tested  and  the  models  merely  swapped 
in  and  out  of  the  simulation.  The  speckle  model  could  also  be  verified  against  some 
lab  tests  with  hardware. 

Real-World  Data.  The  existing  simulation  could  be  run  again  with  some  real- 
world  image  data  and  the  results  compared  to  the  existing  simulation  data.  Real- world 
data  could  be  gathered  from  surfaces  of  varying  degrees  of  optical  roughness  and  the 
performance  of  the  algorithm  examined  for  each  surface. 
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Appendix  A.  Matlab  Code 


r  j  1  his  appendix  contains  all  the  MATLAB®  code  used  in  the  thesis. 


A .  1  Parent  Code 

Listing  A.l:  The  parent  hie  from  which  all  other  functions  are  called. 

(appendixl/ParentCodeVll.m) 


7.7.  Thes i s _s imul at i on_c ode_V 1 1  .m 
7.  Written  by  FLTLT  Brett  Monz 
4  7.  2008 


1 7. 7. 7. 7. 7. 1 7. 7. 7. 7. 1 1 1 1 1 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 1 1 1 1 1 1 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7  7. 


7. 

7. 

7. 

9  7. 
7. 
7. 


7. 

7. 

14  7. 

7. 


7. 

7. 

7. 

19  7. 

7. 

7. 

7. 

7. 

24  7. 

7. 

7. 


7. 

7. 

29  7. 

7. 

7. 

7. 


7. 

34  7. 
7. 
7. 


Version  Development  notes: 

VI  -  Extended  source  modelling 

V2  -  Added  centroiding  algorithm  and  plot  or  noise  vaiance 
V3  -  corrected  error  in  D  aperture  diameter  now  0.07m 
-  added  SRI  eqn  5.38  Roggemann 

V4  -  changed  beta  =1,  s  =  0.5  changed  code  to  only  claculate  .. 
psf 1  once 

V4A-  removed  taking  real  part  of  fft  of  extended  source 
V5  -  added  projection  vector  method  code 
V6  -  added  eqn  5:30  from  Roggemann 

V6aug21  -  corrected  code  such  that  projection  now  performs  ... 
better  than 

centroid  with  extended  source 
V8  -  Thesis  version 

Changed  name  to  Thesis_simulation_code 

VI  -  implemented  centroid  code  as  function  call 

Via  -  implemented  projection  vector  method  as  function  call 
V2  -  implementing  speckle  image  processing  with  centroiding  NO 
turbulence  consideration  implemented  for  speckle 
V3  -  implementing  speckle  image  processing  for  projection  NO 
turbulence  consideration  implemented  for  speckle 
V3A-  Initial  Image  plots  formatted  for  paper/thesis 
V4  -  Turbulence  Included,  now  uses  speckle_gen2 . m  function  and 
save.plot s 

variable  to  allow  sving  to  . eps  or  not 
V4A-  test  to  change  dimension  of  extended  source 

V5  -  added  set  colour  values  to  specify  plot  colours  for  graphs 
V6  -  zero  tilt  phase  screens 

V 7  -  individual  phase  screen  for  each  realisation 
V8  -  use  average  of  100  images  for  reference  for  projection  ... 
method ,  also 

requires  pr o j ect i on_me thod2 . m 

V10  -  changed  speckle  simulation  in  non  turbulent  case 

VII  -  changed  speckle  simulation  in  turbulent  case 


7. 7. 7. 7. 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
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clear  all ; 
close  all ;  clc 

39  check.image  =  true;  "/,  set  if  image  display  req 

7«  che  ck_  image  =  false;  °i  set  if  image  display  not  req 
7«turbulence_included  =  true;  7»  set  to  true  if  turbulence  ... 

considered  and  will  load  phase  screen  and  use  correct  plot  text 
turbulence_included  =  false;  7.  set  to  false  if  turbulence  not  ... 
considered  and 

7«will  use  correct  plot  text  and  not  load  phase  screen 
44  7«save_plots  =  true;  ‘/.uncomment  to  save  plots  to  .  eps  files 

save.plots  =  false;  7«  uncomment  to  not  save  plots  to  .eps  files 
darkgreen  =  [78  136  49]  .  / 2 5 5 ;  7.  set  colour  dark  green 

maroon  =  [153  0  102]  ,/255;  7«  set  colour  maroon 

7o7«  initialising 
49  if  turbulence.included 

load  Simdata_thesis_1000_tilt_rem  7oSimdata_thesis7oSimdata 

end  7oif  turbulence  included  load  individual  otfs  peviously  created 

N  =  128;  7.number  of  aperture  pixels 
54  D  =  .07;  7.  aperture  diameter  in  m 

P  =  64;  ^pixels  in  aperture 
ro  =  0.07;  7.  fried  seeing  parameter 
ua  =  .75;  °i  atmospheric  visibility  factor 

beta  =  .5;  7.  angular  size  of  beacon  relative  to  seeing  limited  ... 

angle 

59  s  =  0.5  ;  7.  normalised  shear 

speckle_source_dim  =  4;  7. dimension  of  source  for  speckle  image 
images  =  1000;  7.  number  of  images  per  light  level 

Kbar  =  (100:100:1000);  7.  set  array  of  kbar  values 
speckle.ref .image  =  zeros  (128); 

64  pro j _ref .image  =  zeros(128); 

phase_slope_centroid_speckle  =  zeros  (1  ,  length  (Kbar ))  ;  7«create  ... 
array  for  phase  slope  values 

phase_slope_centroid_noise  =  zeros  (  1  ,  length  (  Kbar  ))  ;  7«create  array 
for  phase  slope  values 

sigma.n  =  zeros  (1  ,  length  (Kbar)  )  ;  7»create  array  for  tilt  error  ... 
values 

69 

extended.source  =  zeros  (  128 , 128)  ;  7« 
extended.source (64 : 67 , 64 : 67)  =  ones (4 , 4)  ; % 

Likelihood  =  zeros(l,10); 

sigma_n_pro j ection  =  zeros ( 1 , length ( Kbar )) ; 

74  sigma_n_centroid_speckle  =  zeros ( 1 , length ( Kbar )) ; 
sigma_n_pro j ection.speckle  =  zeros ( 1 , length ( Kbar )) ; 
phase_slope_pro j ection  =  zeros  (  1  ,  100) ; 
phase_slope_pro j ection.speckle  =  zeros ( 1  ,  100) ; 


79 

otf_no_turb  =  Make_otf  (32 , 0 , 128 , 1  ,  zeros  ( 128 , 128)  )  ;  7.make  otf 
without  atmospheric  turbulence  screen 
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if  turbulence_included 

otfl  =  otf_array(:  ,  :  ,1)  ;  7»  choose  off  from  tilt  removed  phase 

screen 

84  else 

otfl  =  otf_no_turb; 

end  7.  if  turbulence  included  load  phase  screen  data 

"/«  if  check.image  debug  flag  set  display  images 
89  if  check.image 

psfl  =  ( abs ( if f 1 2 ( ( f f 1 2 ( ext  ended. s our c e ) )  . *  otfl))); 

psf  l_correct_kbar  =  psfl  *  Kbar(5);  "/«  set  average  rxd  photon 
count  Kbar 

temp=ones(128,128) ,/(ones(128,128)+psfl_correct_kbar/2) ; 
speck_image_correct_Kbar  =  icdf(’nbin’,rand(128,128),l, temp ) ; 
i  set  average  rxd  photon  count  Kbar 
94 

7  speck.image  =  speckle_gen2(speckle_source_dim,... 

turbulence_included  ,  otf  1 )  ;  '/.generate  speckle  image 
"/«  7npsfla  =  fftshift(abs(ifft2(otfl))); 

7.  psfl  =  (  abs  (  if  f  1 2  (  (  f  f  1 2  (  ext  ended,  sour  ce  )  )  .  *  otfl))); 

99  set(0,  ’ def aultt ext int erpr et er 1 ,  ’latex’); 
f  =f igure ()  ; 

ax2  =  axe s (’ Unit s ’ ,  ’Inches’,  ’ OuterPosition ’ ,  [015  5]); 

imagesc(psfl (54:77,54:77)) ; colorbar 
colormap( ’gray’) ; 

104  cmap  =  colormap ( ’ gray ’ ) ; 

colormap  (  f  lipud  (  cmap  )  )  ;  7«  use  inverse  colormap  for  printing 

grid 

xlabel([’  X  Dimension  (pixels)  ’  ]); 

ylabel([’Y  Dimension  (pixels)’]); 

109 

f  =  f igure ()  ; 

ax3  =  axe s (’ Unit s ’ ,  ’Inches’,  ’OuterPosition’,  [015  5]); 
imagesc(speck_image_correct_Kbar (54:77,54:77)) ; colorbar 
colormap ( ’gray’) ; 

114  cmap  =  colormap (’ gray ’) ; 

colormap  (  f  1  ipud  (  cmap  ))  ;  7«  use  inverse  colormap  for  printing 

grid 

xlabel([’  X  Dimension  (pixels)  ’  ]); 

ylabel([’Y  Dimension  (pixels)’]); 

119  end  7,  if  check.image 

7.7.  Make  Images  at  "images"  light  levels 

for  iloop2  =  1  :  length  ( Kbar )  7.  loop  for  light  levels 

124 

for  iloopl  =  1:  images  7o  loop  through  images  per  light  level 
iloopl 
pause ( . 1 ) 

if  turbulence.included 
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129  otfl  =  otf  .array  (  :  ,  :  ,  iloopl )  ;  7.  choose  otf  from  tilt 

removed  phase  screen 
if  iloopl  ==  1 

for  iloop3  =  1: images 

psfl  =  ( abs ( if f 1 2 ( ( f f 1 2 ( ext  ended. sour ce ) )  . * 
otfl)))  ; 

psf  l_correct_kbar  =  psfl  *  Kbar  ( iloop2 )  ;  7o  set 
average  rxd  photon  count  Kbar 

134  temp  =  ones (128 , 128)  ./( ones ( 128 , 128) +  ..  . 

psf l_correct_kbar/2) ; 

speck_image_correct_Kbar  =  icdf ( ’ nbin ’  , rand  .  .  . 
(128,128)  ,1,  temp)  ;  7»  set  average  rxd  photon 

count  Kbar 

speckle.ref .image  =  speckle.ref .image  +  (... 

speck_image_correct_Kbar) /images ; 
pro j _ref .image  =  pr oj _ref _ image  +  (... 
psf l_correct_ kbar)/ images ; 

139  end  7,  for  iloop3 

end  7.  if  iloopl 

else 

otfl  =  otf_no_turb; 

144  pro j _ref _image=psf 1 *Kbar ( iloop2 ) ; 

speckle_ref_image  =  psfl*Kbar(iloop2)  ; 
end  7oif  turbulence  included  load  phase  screen  data 
7psfla  =  fftshift(abs(ifft2(otfl))); 

psfl  =  ( abs ( if f 1 2 ( ( f f 1 2 ( ext  ended. sour ce ) )  . *  otfl))); 

149 

psf  l_correct_kbar  =  psfl  *  Kbar  ( iloop2 )  ;  °L  set  average  rxd 
photon  count  Kbar 

t emp  =  one  s (128,128)  . /(ones (128, 128) +psfl_correct_kbar/2)  ; 
speck_image_correct_Kbar  =  i cdf ( ’ nbin’, rand(128, 128),  1,... 

temp);  7»  set  average  rxd  photon  count  Kbar 
psf _noise_less  =  psf l_correct_kbar ; 

154  Cx_no_noise  =  cent r o id ( ext ended.sour ce ) ; %  Cx_no_noise  =  ... 

cent r o id ( psf 1 _cor r ect _kbar ) ; 

Cx_no_noise_speckle  =  c  ent  r  o  id  (  ext  ended,  sour  ce  )  ;  7o  ... 

Cx_no_noise_speckle  =  centroid( speck_image_correct_Kbar ) 

i 

psf  l_correct_kbar_noisy  =  po  i  s  srnd  (  psf  1  _c  or  r  e  ct  _kbar  )  ;  7.  . 

make  noisy 

t emp  =  one  s (128,128)  ,/(ones(128,128)+psfl_correct_kbar)  ; 

159  speck_image_correct_Kbar_noisy=  icdf ( 1 nbin 1 , rand ( 128 , 128) . 

,  l.temp);  ;  7«  make  noisy 

7,  calculate  centroid  noisy  for  extend  source 
Cx.noise  =  cent r o id ( psf 1 _ c or r e ct _kbar _no i sy ) ; 
phase_slope_centroid_noise (iloopl)=  (Cx_noise-Cx_no_noise) 
*2*pi*P/(N*D);  7.  in  rads/m 

164 
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%  calculate  centroid  noisy  for  speckle  source 
Cx_no_noise_speckle  =  c ent r o id ( spe ck_ image. corr e ct _Kbar ) ; 
Cx_noise_speckle  =  centroidC speck_image_correct_Kbar_noisy 
)  ; 

169  phase_slope_centroid_speckle ( iloopl ) =  (Cx_noise_speckle  - 

Cx_no_noise)  *2*pi*P/ (N*D)  ;  7.  in  rads/m 

7,  calculate  phase  slope  using  projection  method 
phase_slope_projection( iloopl )=  projection_method3(. .  . 

psf  1  _  corre  ct  _kbar_no  i  sy  ,  pro  j  _ref  _  image  ,  P  ,  N  ,  D)  ;  7»  in  rads 
/  m 

174  phase_slope_proj ection_speckle ( iloopl ) =  proj ection_method2 

(speck_image_correct_Kbar,speckle_ref_image,P,N,D)  ;7«  in 
r ads /m 

end  7»  iloopl 

s igma_n ( iloop2 )  =  sqrt ( sum ( ( phase. slope. cent r o id_no i se ) . ~ 2 ) . / . 
images ) ; 

s igma_n_pr o j e ct i on ( i loop2 )  =  sqrt ( sum ( (phase_slope_proj ection) 

. ~2) .  / images ) ; 

179  sigma_n_centroid_speckle ( iloop2 )  =  sqrt ( sum ((..  . 

phase_slope_centroid_speckle) .  ~2) ./images) ; 

sigma_n_pro j ection_speckle ( iloop2 )  =  sqrt ( sum ((..  . 
phase_slope_projection_speckle) . ~ 2 ) ./images) ; 


iloop2 

184  end  7«  iloop2 
7.7.  plot  figure 

set(0,  ’ def aultt ext int erpr et er  ’  ,  ’latex'); 
f  =f igure ()  ; 

189  axl  =  axe s (’ Unit s ’ ,  ’Inches’,  ’ OuterPosition ’ ,  [015  5]); 

plot ( axl  , Kbar , sigma.n ,  ’ b- -o  ’  )  ; 
hold  on;  grid  on; 

plot ( axl  , Kbar , s igma_n_pr o j ection,  ’ r - . * ’ )  ; 

194  plot ( axl  , Kbar , sigma_n_ cent roid_ speckle  ,  ’ - . d ’  ,  ’color’  ,darkgreen)  ; 
plot ( axl  , Kbar , s igma_n_pr o j  ection.speckle ,  ’ - . p ’  ,  ’color’  , maroon)  ; 
xlabel([’  $\bar{K}  $  -  Average  photons  received  per  image  ’  ]); 

ylabel ( [ ’ Tilt  Error  Due  to  Noise  (rad/m)’]); 

199  12  =  legend([’Centroiding  method’], [’Projection  vector  method’], [’ 
Centroid  method  on  Speckle ’],[’ Proj ection  method  on  Speckle’]); 
set(12,  ’location’,  ’northeast’,  ’Fontsize’, 10  ); 

if  save.plots 

if  turbulence.included 

204  print (’ -depsc ’  ,  ’ Research_tilt_error_turbulence  ’  )  ; 

else 
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209 


print ( ’ -depsc ’  ,  ’ Research._tilt_error_no_turbulence  1  )  ; 

end  °/0  if  turbulence. included) 
end  y.  if  save.plots 

'/.save  (  ’  Thesis_sim_data_no_turbulence_1000_4x4_ext_source_testingl  .  .  . 
’  ,  ’Kbar’,  1 sigma.n  ’  ,  ’ sigma_n_proj ection  ’  ,  ’  . .  . 
sigma.n.centroid. speckle  1  ,  ’ sigma.n.pro j ection.speckle ’ ) 
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A. 2  Make  Otf  function 


Listing  A. 2:  This  function  computes  the  Optical  transfer  Function  from  the  source 
to  the  lens 

(appendixl /Makeotf2.m) 

function  [otf , aperture]  =  Make_ot f 2 ( rl  ,  r2 , si , scale  ,  phase ) 

%  [otf , apeture]  =  make.otf (rl ,r2 , si , scale , phase) ; 

mi  =  floor(si/2); 

5  mi =mi + 1 ; 

aperture  =  zeros ( si , si )  ; 
for  i  =  1 : si 

for  j  =  1 : si 

10 

dist  =  sqrt ( (i-mi) ~2+ ( j -mi) ~2) ; 
if (dist  <  =  rl ) 

if (dist  >  =  r2) 

15  apertur e ( i  ,  j )  =  1; 


end 

end 


20  end 


end 

H  tilt  removal  section 
tx  =  2*pi*(-64:63) / 128 ; 

25  tx.matrix  =  one s ( 1 28 , 1 ) * tx ; 

wvs_x  =  sum ( sum ( tx.matr ix . *  phase .* apertur e ))  ; 
wvs_x  =  wvs_x/ sum ( sum ( apertur e .* tx.matrix . ~ 2 )) ; 
new_screen  =  phase  -  wvs_x*tx_matrix ; 

V/ 

30  pupil  =  apert ur e . *  cos ( new_s cr een)  +  sqrt (- 1 )* apertur e .* s in (..  . 
new_screen) ; 

psf  =  r eal ( f f 1 2 ( pupi 1 )  . *  con j  ( f f 1 2 ( pupi 1 ) ) )  ; 
psf  =  s cale *psf / sum ( sum ( psf ) ) ; 

35  otf  =  fft2(psf); 
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A. 3  Centroiding  Algorithm 


Listing  A. 3:  This  function  computes  the  centroid  of  the  image  sent  and  returns  the 
wavefront  tilt. 

(appendixl  /  centroid,  m) 

funct ion[cx]=centroid(image) 

%  returns  centroid  (cx)  in  x  dimension  only  of  (image) 

V.  image  should  have  dim  128x128  (not  checked) 

[posx,posy]  =  me  shgr  id  (  -  1 2  :  1  :  1 1 )  ;  °/«  create  position  array  each  ... 

element  value  is  its  delta  from  center  in  pixels 
5  cx=  sum ( sum (posx.*image (54:77,54:77)))  ./ sum ( sum (image (54:77,54:77))  .  .  . 

)  ; 
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A. 4  Projection  Method  Algorithm  2 


Listing  A. 4:  This  function  computes  the  wavefront  tilt  using  a  MLE  algorithm  with 

fractional  and  integer  pixel  resoultion. 

(appendixl/projectionmethod2.m) 

funct ion[phase_slope]=projection_method2(image ,ref_image ,P,N,D) 

1  returns  phase  slope  of  (image) 

7.  N  number  of  aperture  pixels 
°i  D  aperture  diameter  in  m 
5  */.  P  pixels  in  aperture 

•/.  image  should  have  dim  128x128  (not  checked) 

7,  requires  makeshift_vec.m  to  reside  in  same  directory 
bvec  = -5  :  .  01  :  5  ;  7»  vector  of  pixel  locations  for  linear  ... 
interpolation 

psf 11=( image (54:77, 54:77)/2) ; 

10  Px  =  sum(psfll); 

Po  =  sum ( r ef _ image ( 54 : 77 , 54  :  77 ) /2 )  ; 
bindex  =  1  ; 

7,  do  linear  interpolation 
for  bloop  =  -5:. 01:5 
15  bi  =  floor (bloop) ; 

bf  =  bloop -bi  ; 

Pso  =  makeshif t_vec (Po , bi )  ; 

Psl  =  makeshif t_vec (Po , bi+1) ; 

Pos  =  (l-bf)*Pso  +  bf*Psl; 

20 

7.  do  likelihood 

Likelihood (bindex)  =  sum ( Px . *  log ( Pos ) -  Pos); 
bindex  =  bindex  +  1  ;  7«  update  bindex 

end7.bloop 

25  maxLikelihood  =  max (Likelihood) ; 

dummyval  =  f ind ( L ike  1 iho od  ==  maxLikelihood); 
phase_slope=  bvec  ( dummyval  )*2*pi*P/(N*D)  ;7»  in  rads/m 
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A. 5  Projection  Method  Algorithm  3 

Listing  A. 5:  This  function  computes  the  wavefront  tilt  using  a  MLE  algorithm  . 

(appendixl/projectionmethod3.m) 

funct ion[phase_slope]=projection_method3(image , ref .image ,P,N,D) 

7.  returns  phase  slope  of  (image) 

3  7„  N  number  of  aperture  pixels 
7.  D  aperture  diameter  in  m 
7.  P  pixels  in  aperture 

7.  image  should  have  dim  128x128  (not  checked) 

"/«  requires  makeshift_vec.m  to  reside  in  same  directory 
8  bvec  = -5  :  .  01  :  5  ;  7.  vector  of  pixel  locations  for  linear  ... 

interpolation 

psf 11=( image (54:77, 54:77)/2) ; 

Px  =  sum ( psf 1 1 )  ; 

Po  =  sum ( r ef _ image ( 54 : 77 , 54 : 77 ) /2 ) ; 
bindex  =  1  ; 

13  7.  do  linear  interpolation 
for  bloop  =  -5:. 01:5 
‘/,  bi  =  floor  (bloop)  ; 

7.  bf  =  bloop-bi; 

7.  Pso  =  makeshif  t_vec  (Po  ,  bi)  ; 

18  7.  Psl  =  make shif  t _ ve c  ( Po  ,  bi  + 1 )  ; 

Pos  =  makeshif t_vec (Po , bloop ) ; 

7.  do  likelihood 

Likelihood (bindex)  =  sum ( Px . *  log ( Pos ) -  Pos); 

23  bindex  =  bindex  +  1  ;  7«  update  bindex 

end7.bloop 

maxLikelihood  =  max (Likelihood) ; 

dummyval  =  f ind ( L ike  1 ihood  ==  maxLikelihood); 

phase_slope=  bvec ( dummyval )*2*pi*P/(N*D)  ; %  in  rads/m 
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A. 6  Vector  Linear  Shift  Algorithm 

Listing  A. 6:  This  function  implements  a  linear  shift  of  points  in  a  supplied  vector. 

(appendixl  /  makeshift  vec.rn) 

function  [img2]=  makeshif t_vec ( imgl , dx) 

%  function  [img2]=  makeshif t ( imgl , dx , dy) 

3  */.  dy  and  dx  are  the  shifts  in  the  vertical  and  horizontal  ... 
directions  respectively 

V.  imgl  and  img2  are  the  two  images  from  a  sequence  of  video 

'/.  delta  is  the  denominator  of  the  fraction  of  a  pixel  to  which  the... 

estimation  is  to  be  done 
*/.  ex  1/10  pixel  estimation  means  delta  =10 

8  sz= s ize ( imgl ) ; 
sz=max(max(sz) ) ; 
center  =  [f loor ( sz/ (2) ) +1] ; 

linx  =  -  cent  er  +  1  :  1  :  cent  er -2  ; 
linx  =  -2* pi  *  1 inx/ sz  ; 

13  linx  =  f f t shi f t ( 1 inx ) ; 


px  =  cos ( 1 inx*dx ) + sqrt ( - 1 ) * s in ( 1 inx * dx ) ; 
18  img2  =  real ( if f t ( f f t ( imgl )  . * (px) ) )  ; 
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A.  7  Atmospheric  Turbulence  Simulation 

Listing  A. 7:  This  code  produces  tilt  removed  phase  screen  realizations. 

(appendixl  /PhaseScreen.m) 

1  function  [screen , deltrho] =PHASE_SCREEN (D , ro , laml , si , K) 

7,  [screen]  =  PHASE_SCREEN  (D  ,  ro  ,si,K); 

°/oD  is  the  linear  dimension  in  meters  corresponding  to  a  square  ... 
with  si  pixels  on  a 

‘/.side .  Example:  You  have  a  wfs  aperture  of  10  cm  and  are  putting 

it 

1  in  an  array  32  pixels  on  a  size.  The  array  is  twice  the  size  of 
the  aperture . 

6  */.  then  D=.l  and  si  =  32.  ro  is  Fried’s  seeing  parameter  in  meters 

1  K  is  the  number  of  independent  phase  screens  you  need  to  ... 
generate . 

°/0on  output  make  sure  deltrho  is  not  larger  than  about  4  ... 
millimeters 

s i 1 = 1 024 ; 

11  si2=1024 ; 

deltrho  =  D/ si  ; 

Lo=deltrho*sil/2 
ac  =  zeros (sil , si2)  ; 
mi 1  =  si 1 /2+ 1 ; 

16  mi2=si2/2+l ; 
for  i  =  1 : si 1 

i  ; 

for  j  =  l:si2 

rho  =  delt rho* sqrt ( ( i -mi 1 ) “ 2+ ( j -mi  2 ) “ 2)  ; 

21  ac (i , j ) =besselk ( (5/6) , (2*pi* (rho) /Lo) ) * ( (rho) “ (5/6) ) * (Lo . . 

/ (2*pi) ) “ (5/6)  ; 

ac  (  i  ,  j  )  =ac(i,j)/ ( (2~ (5/6) ) *  gamma  (11/6)  )  ; 
if (rho==0) 

ac  (i  ,  j  )  =0  ; 

end 

26  end 

end 

sz  =  1 024 ; 

Cn2dz  =  ( ( . 185) “ (5/3) ) *(laml“2) /(ro“ (5/3) ) *10“ (-12)  ; 
k  =  2*pi/(laml*10“(-6)); 

31  acl=ac* . 033* (4*pi“2) *k*k*Cn2dz ; 
acl=fftshift(acl) ; 
ac 1 (1 , 1 ) =  ac 1  (1,2); 
ftl  =  sqrt ( abs ( f f t 2 ( ac 1 ) ) ) ; 
for  k= 1 : K 

36  phase  =  r andn ( sz , sz )  ; 

phasel  =  r eal ( if f 1 2 ( f f 1 2 ( phase )  .  * f 1 1 )  )  ; 
screen(: ,: ,k)=phasel(l:si,l:si); 
k 

pause ( . 1 ) 

41  end ; return 
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word.  Page  numbers  in  bold  represent  concept  definition  or  introduction. 


AO,  2 
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infra  red,  see  IR 

IR,  1 

maximum  likelihood  estimator,  see  MLE 

ML,  30 

MLE,  1 

01,  2 

OPL,  26 

optical  path  length,  see  OPL 

OTF,  38 

PDF,  17 

Point  Spread  Function,  see  PSF 

PSF,  6,  38 

RV,  18 
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